xmdsintegratesiip.cc

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

  // 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,"");

  list<XMLString> myNonLoopPropagationCodeList=*nonLoopPropagationCodeList();   
  list<integrateMomentGroup> myIntegrateMomentGroupList=*integrateMomentGroupList();
    
  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 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");

      // 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,"\n");


      char indent[64];
      for(long unsigned int i=0;i<nDims;i++) {
	indent[i]=0x09;
      }
      indent[nDims]=0;

      if((simulation()->parameters()->stochastic)&&(!noNoises())) {
	if(simulation()->parameters()->errorCheck) {
	  fprintf(outfile,"%sif(_generator_flag==1)\n",indent);
	  fprintf(outfile,"%s	_make_noises(_gen1,_var,_noises,_n_noises);\n",indent);
	  fprintf(outfile,"%selse if(_generator_flag==2)\n",indent);
	  fprintf(outfile,"%s	_make_noises(_gen2,_var,_noises,_n_noises);\n",indent);
	  fprintf(outfile,"%selse {\n",indent);
	  fprintf(outfile,"\n");
	  fprintf(outfile,"%s	_make_noises(_gen1,_var/2,_noises,_n_noises);\n",indent);
	  fprintf(outfile,"%s	_make_noises(_gen2,_var/2,_noises2,_n_noises);\n",indent);
	  fprintf(outfile,"\n");
	  fprintf(outfile,"%s	for(unsigned long _s0=0;_s0<_n_noises;_s0++)\n",indent);
	  fprintf(outfile,"%s		_noises[_s0] += _noises2[_s0];\n",indent);
	  fprintf(outfile,"%s	}\n",indent);
	  fprintf(outfile,"\n");
	}
	else {
	  fprintf(outfile,"%s_make_noises(_gen,_var,_noises,_n_noises);\n",indent);
	  fprintf(outfile,"\n");
	}
      }

      fprintf(outfile,"%sfor(unsigned long _i%li=0;_i%li<_%s_main_ncomponents;_i%li++)\n",indent,nDims+1,nDims+1,fieldName,nDims+1);
      fprintf(outfile,"%s	_%s_main_old[_i%li] = _active_%s_main[_%s_main_index_pointer + _i%li];\n",indent,fieldName,nDims+1,fieldName,fieldName,nDims+1);
      fprintf(outfile,"\n");

      if(crossVectorList.size() > 0) {

	for(list<const xmdsVector*>::const_iterator ppxmdsVector = crossVectorList.begin(); ppxmdsVector != crossVectorList.end(); ppxmdsVector++) {

	  const char* cvName = (*ppxmdsVector)->name()->c_str();

	  fprintf(outfile,"%sfor(unsigned long _i%li=0;_i%li<_%s_%s_ncomponents;_i%li++)\n",
		  indent,nDims+1,nDims+1,fieldName,cvName,nDims+1);
	  fprintf(outfile,"%s	_%s_%s_old[_i%li] = _active_%s_%s[_%s_%s_index_pointer + _i%li];\n" ,
		  indent,fieldName,cvName,nDims+1,fieldName,cvName,fieldName,cvName,nDims+1);
	  fprintf(outfile,"\n");
	}
      }

      if(nIterations()>1) {
	fprintf(outfile,"%sfor(unsigned long _i%li=0;_i%li<%li;_i%li++) {\n",indent,nDims+1,nDims+1,nIterations()-1,nDims+1);
	fprintf(outfile,"\n");
	fprintf(outfile,"// *** propagation (and cross_propagation code) ***\n");
	fprintf(outfile,"%s\n",propagationCode()->c_str());
	if(crossVectorList.size() > 0) {
	  fprintf(outfile,"%s\n",crossPropagationCode()->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] = _%s_main_old[%li] + d%s_d%s*(_step/2);\n",
		  indent,fieldName,fieldName,i,fieldName,i,
		  mainVector->componentName(i)->c_str(),simulation()->parameters()->propDimName.c_str());
	}

	for(list<const xmdsVector*>::const_iterator ppxmdsVector = crossVectorList.begin(); ppxmdsVector != crossVectorList.end(); ppxmdsVector++) {

	  const char* cvName = (*ppxmdsVector)->name()->c_str();

	  fprintf(outfile,"\n");
	  for(long unsigned int i=0;i<(*ppxmdsVector)->nComponents();i++) {
	    fprintf(outfile,"%s	_active_%s_%s[_%s_%s_index_pointer + %li] = _%s_%s_old[%li] + d%s_d%s*(_%s_dx%li/2);\n",
		    indent,fieldName,cvName,fieldName,cvName,i,fieldName,cvName,i,(*ppxmdsVector)->componentName(i)->c_str(),
		    simulation()->field()->geometry()->dimension(crossDimNumber())->name.c_str(),fieldName,crossDimNumber());
	  }
	}
	fprintf(outfile,"%s	}\n",indent);
	fprintf(outfile,"\n");
      }

      fprintf(outfile,"// *** propagation (and cross_propagation code) ***\n");
      fprintf(outfile,"%s\n",propagationCode()->c_str());
      if(crossVectorList.size() > 0) {
	fprintf(outfile,"%s\n",crossPropagationCode()->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] = _%s_main_old[%li] + d%s_d%s*_step;\n",
		indent,fieldName,fieldName,i,fieldName,i,
		mainVector->componentName(i)->c_str(),simulation()->parameters()->propDimName.c_str());
      }

      for(list<const xmdsVector*>::const_iterator ppxmdsVector = crossVectorList.begin(); ppxmdsVector != crossVectorList.end(); ppxmdsVector++) {

	const char* cvName = (*ppxmdsVector)->name()->c_str();

	fprintf(outfile,"\n");
	for(long unsigned int i=0;i<(*ppxmdsVector)->nComponents();i++) {
	  fprintf(outfile,"%s_active_%s_%s[_%s_%s_index_pointer + %li] = _%s_%s_old[%li] + d%s_d%s*(_%s_dx%li/2);\n",
		  indent,fieldName,cvName,fieldName,cvName,i,fieldName,cvName,i,(*ppxmdsVector)->componentName(i)->c_str(),
		  simulation()->field()->geometry()->dimension(crossDimNumber())->name.c_str(),fieldName,crossDimNumber());
	}
      }

      if(crossVectorList.size() > 0) {
	fprintf(outfile,"\n");
	fprintf(outfile,"%sif(_i%li<_main_lattice%li-1) {\n",indent,crossDimNumber(),crossDimNumber());

	for(list<const xmdsVector*>::const_iterator ppxmdsVector = crossVectorList.begin(); ppxmdsVector != crossVectorList.end(); ppxmdsVector++) {

	  const char* cvName = (*ppxmdsVector)->name()->c_str();

	  fprintf(outfile,"\n");
	  fprintf(outfile,"%s	for(unsigned long _i%li=0;_i%li<_%s_%s_ncomponents;_i%li++)\n",
		  indent,nDims+1,nDims+1,fieldName,cvName,nDims+1);
	  fprintf(outfile,"%s		_active_%s_%s[_%s_%s_index_pointer + _i%li + _%s_%s_ncomponents",
		  indent,fieldName,cvName,fieldName,cvName,nDims+1,fieldName,cvName);
	  for(long unsigned int i=crossDimNumber()+1;i<nDims;i++) {
	    fprintf(outfile,"*_%s_lattice%li",fieldName,i);
	  }
	  fprintf(outfile,"] = 2*_active_%s_%s[_%s_%s_index_pointer + _i%li] - _%s_%s_old[_i%li];\n",
		  fieldName,cvName,fieldName,cvName,nDims+1,fieldName,cvName,nDims+1);
	}
	fprintf(outfile,"%s	}\n",indent);
      }

      simulation()->field()->closeLoops(outfile,0,myTotalVectorsList);

    }
    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++;