xmdsmomentgroup.cc

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

	  fprintf(outfile,"	");
	  fprintf(outfile,"if(k%s>((_%s_lattice%li-1)/2 + 0.1)*_%s_dk1) {\n",simulation()->field()->geometry()->dimension(1)->name.c_str(),mgFieldName,main2PostDim(1),fieldName);
	  for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
	    fprintf(outfile,"	");
	    fprintf(outfile,"	_%s_%s_index_pointer += (_%s_lattice1-_%s_lattice%li)",
		    fieldName,pXMLString->c_str(),fieldName,mgFieldName,main2PostDim(1));
	    fprintf(outfile,"*_%s_lattice0",fieldName);
	    for(j=2;j<nDims;j++) {
	      fprintf(outfile,"*_%s_lattice%li",fieldName,j);
	    }
	    fprintf(outfile,"*_%s_%s_ncomponents;\n",fieldName,pXMLString->c_str());
	  }
	  fprintf(outfile,"	");
	  fprintf(outfile,"	k%s -= _%s_lattice%li*_%s_dk1;\n",simulation()->field()->geometry()->dimension(1)->name.c_str(),mgFieldName,main2PostDim(1),fieldName);
	  fprintf(outfile,"	");
	  fprintf(outfile,"	}\n");
	}

        // even if cross-section
        fprintf(outfile,"	");
        fprintf(outfile,"}\n");
    }
    else {
        for(i=nDims;i>0;i--) {
        
	  const char* dimName = simulation()->field()->geometry()->dimension(i-1)->name.c_str();
        
	  fprintf(outfile,"\n");
	  if(*psamplingLatticeIr==0) {
	    // integration
	    if(i==nDims) {
	      for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
		for(j=0;j<i;j++) {
		    fprintf(outfile,"	");
		}
		fprintf(outfile,"_%s_%s_index_pointer += _%s_%s_ncomponents;\n",
			fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str());
	      }
	      fprintf(outfile,"\n");
	    }
	    if(samplingSpace(i-1)) {
	      for(unsigned long j=0; j<i; j++) {
		  fprintf(outfile,"	");
	      }
	      fprintf(outfile,"k%s +=	_%s_dk%li;\n",dimName,fieldName,i-1);
	      for(unsigned long j=0; j<i; j++) {
		  fprintf(outfile,"	");
	      }
	      fprintf(outfile,"if(k%s>((_%s_lattice%li-1)/2 + 0.1)*_%s_dk%li)\n",dimName,fieldName,i-1,fieldName,i-1);	
	      for(unsigned long j=0; j<i; j++) {
		  fprintf(outfile,"	");
	      }
	      fprintf(outfile,"	k%s -= _%s_lattice%li*_%s_dk%li;\n",dimName,fieldName,i-1,fieldName,i-1);
	    }
	    else {
	      for(unsigned long j=0; j<i; j++) {
		  fprintf(outfile,"	");
	      }
	      fprintf(outfile,"%s += _%s_dx%li;\n",dimName,fieldName,i-1);
	    }
	  }
	  else if(*psamplingLatticeIr==1) {
	    // cross-section
	    for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
	      for(j=0;j<i;j++) {
		  fprintf(outfile,"	");
	      }
	      fprintf(outfile,"_%s_%s_index_pointer += (_%s_lattice%li",fieldName,pXMLString->c_str(),fieldName,i-1);
	      if(!samplingSpace(i-1)) {
		fprintf(outfile,"-_i%li",i-1);
	      }
	      if(!(i==nDims)) {
		fprintf(outfile,"-1");
	      }
	      fprintf(outfile,")");
	      for(j=i;j<nDims;j++) {
		fprintf(outfile,"*_%s_lattice%li",fieldName,j);
	      }
	      fprintf(outfile,"*_%s_%s_ncomponents;\n",fieldName,pXMLString->c_str());
	    }
	  }
	  else {
	    // normal sampling
	    if(samplingSpace(i-1)) {
	      // narrow k-space window
	      if(i==nDims) {
		for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
		  for(j=0;j<i;j++) {
		      fprintf(outfile,"	");
		  }
		  fprintf(outfile,"_%s_%s_index_pointer += _%s_%s_ncomponents;\n",
			  fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str());
		}
		fprintf(outfile,"\n");
	      }
	      for(j=0; j<i; j++) {
		  fprintf(outfile,"	");
	      }
	      fprintf(outfile,"k%s +=	_%s_dk%li;\n",dimName,fieldName,i-1);
	      for(j=0; j<i; j++) {
		  fprintf(outfile,"	");
	      }
	      fprintf(outfile,"if(k%s>((_%s_lattice%li-1)/2 + 0.1)*_%s_dk%li) {\n",
		      dimName,mgFieldName,main2PostDim(i-1),fieldName,i-1);
	      for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
		for(j=0;j<i;j++) {
		    fprintf(outfile,"	");
		}
		fprintf(outfile,"	_%s_%s_index_pointer += (_%s_lattice%li-_%s_lattice%li)",
			fieldName,pXMLString->c_str(),fieldName,i-1,mgFieldName,main2PostDim(i-1));
		for(j=i;j<nDims;j++) {
		  fprintf(outfile,"*_%s_lattice%li",fieldName,j);
		}
		fprintf(outfile,"*_%s_%s_ncomponents;\n",fieldName,pXMLString->c_str());
	      }
	      for(j=0; j<i; j++) {
		  fprintf(outfile,"	");
	      }
	      fprintf(outfile,"	k%s -= _%s_lattice%li*_%s_dk%li;\n",dimName,mgFieldName,main2PostDim(i-1),fieldName,i-1);
	      for(j=0;j<i;j++) {
		  fprintf(outfile,"	");
	      }
	      fprintf(outfile,"	}\n");
	    }
	    else {
	      // coarse x-space sampling
	      for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
		for(j=0;j<i;j++) {
		    fprintf(outfile,"	");
		}
		fprintf(outfile,"_%s_%s_index_pointer += (_%s_lattice%li/_%s_lattice%li",
			fieldName,pXMLString->c_str(),fieldName,i-1,mgFieldName,main2PostDim(i-1));
		if(!(i==nDims)) {
		  fprintf(outfile,"-1");
		}
		fprintf(outfile,")");
		for(j=i;j<nDims;j++) {
		  fprintf(outfile,"*_%s_lattice%li",fieldName,j);
		}
		fprintf(outfile,"*_%s_%s_ncomponents;\n",fieldName,pXMLString->c_str());
	      }
	      fprintf(outfile,"\n");
	      for(j=0;j<i;j++) {
		  fprintf(outfile,"	");
	      }
	      fprintf(outfile,"%s += _%s_dx%li;\n",dimName,mgFieldName,main2PostDim(i-1));
	    }
	  }
        
	  if(i==1&simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
	    // even if cross-section
	    fprintf(outfile,"	");
	    fprintf(outfile,"}\n");
	  }
	  else
	    if(!(*psamplingLatticeIr==1)) {
	      // if not cross-section
	      for(j=0;j<i;j++) {
		  fprintf(outfile,"	");
	      }
	      fprintf(outfile,"}\n");
	    }
        
	  psamplingLatticeIr++;
	}
    }
  }

  fprintf(outfile,"\n");

  fprintf(outfile,"_%s_x0[_%s_sample_pointer]=%s;\n",mgFieldName,mgFieldName,simulation()->parameters()->propDimName.c_str());
  fprintf(outfile,"\n");

  if(simulation()->parameters()->nPaths==1) {
    if(simulation()->parameters()->usempi) {
      fprintf(outfile,"printf(\"Rank[%%i]  Sampled field (for moment group #%li) at %s	= %%e\\n\",rank,%s);\n",
	      myGroupNumber+1,simulation()->parameters()->propDimName.c_str(),simulation()->parameters()->propDimName.c_str());
    }
    else {
      fprintf(outfile,"printf(\"Sampled field (for moment group #%li) at %s	= %%e\\n\",%s);\n",
	      myGroupNumber+1,simulation()->parameters()->propDimName.c_str(),simulation()->parameters()->propDimName.c_str());
    }
    fprintf(outfile,"\n");
  }
  fprintf(outfile,"_%s_sample_pointer++;\n",mgFieldName);
  fprintf(outfile,"}\n");
  fprintf(outfile,"\n");

  // ***************************************
  // ******** process routine  *************
  // ***************************************

  fprintf(outfile,"// *************************\n");
  fprintf(outfile,"void _%s_process() {\n",mgFieldName);
  fprintf(outfile,"\n");
  fprintf(outfile,"double *_to_field;\n");
  if(simulation()->parameters()->nPaths>1) {
    fprintf(outfile,"double *_to_field_sd;\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());
    }
  }

  if(myPostCode.length()>0) {
    for(list<XMLString>::const_iterator pXMLString = myPostMomentsList.begin(); pXMLString != myPostMomentsList.end(); pXMLString++) {
      fprintf(outfile,"double %s;\n",pXMLString->c_str());
    }
  }

  if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
      if(complexRawVector) {
        fprintf(outfile,"complex* _temp_vector;\n");
      }
      else {
        fprintf(outfile,"double* _temp_vector;\n");
      }
      fprintf(outfile,"unsigned long _length;\n\n");

      fprintf(outfile,"_length = _%s_size*_%s_fullstep_ncomponents;\n",mgFieldName,mgFieldName);
      if(complexRawVector) {
        fprintf(outfile,"_temp_vector = new complex[_length];\n");
      }
      else {
        fprintf(outfile,"_temp_vector = new double[_length];\n");
      }
    
      if(complexRawVector) {
        fprintf(outfile,"MPI_Reduce(_%s_raw,_temp_vector,2*_length,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);\n\n",mgFieldName);
      }
      else {
        fprintf(outfile,"MPI_Reduce(_%s_raw,_temp_vector,_length,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);\n\n",mgFieldName);
      }
    
      fprintf(outfile,"for(unsigned long _i0=0; _i0<_length; _i0++)\n");
      fprintf(outfile,"        _%s_raw[_i0]=_temp_vector[_i0];\n",mgFieldName);
      fprintf(outfile,"\n");
      fprintf(outfile,"delete[] _temp_vector;\n\n");
    }

  fprintf(outfile,"\n");

  /* if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic)
     {
     fprintf(outfile,"if(rank==0) {\n");
     } */

  // assign pointers
  if(simulation()->parameters()->errorCheck) {
    fprintf(outfile,"if(_half_step) {\n");
    fprintf(outfile,"	_to_field = _%s_halfstep;\n",mgFieldName);
    if(simulation()->parameters()->nPaths>1) {
      fprintf(outfile,"	_to_field_sd = _%s_halfstep_sd;\n",mgFieldName);
    }
    fprintf(outfile,"	}\n");
    fprintf(outfile,"else {\n");
    fprintf(outfile,"	_to_field = _%s_fullstep;\n",mgFieldName);
    if(simulation()->parameters()->nPaths>1) {
      fprintf(outfile,"	_to_field_sd = _%s_fullstep_sd;\n",mgFieldName);
    }
    fprintf(outfile,"	}\n");
  }
  else {
    fprintf(outfile,"_to_field = _%s_fullstep;\n",mgFieldName);
    if(simulation()->parameters()->nPaths>1) {
      fprintf(outfile,"_to_field_sd = _%s_fullstep_sd;\n",mgFieldName);
    }
  }
  fprintf(outfile,"\n");

  vectors2space(outfile,myPostSpace,myPostProcessingVectorNamesList,"");

  fprintf(outfile,"unsigned long _%s_raw_index_pointer=0;\n",mgFieldName);
  fprintf(outfile,"unsigned long _%s_processed_index_pointer=0;\n",mgFieldName);
  fprintf(outfile,"\n");

  if(postSpace(0)) {
    fprintf(outfile,"double k%s = 0;\n",geometry()->dimension(0)->name.c_str());
  }

  fprintf(outfile,
	  "for(unsigned long _i0=0; _i0<_%s_lattice0; _i0++) {\n",mgFieldName);
  fprintf(outfile,"\n");

  if(!postSpace(0)) {
    fprintf(outfile,"%s = _%s_x0[_i0];\n",geometry()->dimension(0)->name.c_str(),mgFieldName);
  }
  fprintf(outfile,"\n");

  for(i=1; i<myNDims; i++) {
    for(unsigned long j=0; j<i; j++) {
	fprintf(outfile,"	"); 
    }
    if(postSpace(i)) {
      fprintf(outfile,"double k%s = 0;\n",geometry()->dimension(i)->name.c_str());
    }
    else {
      fprintf(outfile,"double %s = _%s_xmin%li;\n",geometry()->dimension(i)->name.c_str(),mgFieldName,i);
    }
    fprintf(outfile,"\n");
    for(unsigned long j=0; j<i; j++) {
	fprintf(outfile,"	");
    }
    fprintf(outfile,"for(unsigned long _i%li=0; _i%li<_%s_lattice%li; _i%li++) {\n",i,i,mgFieldName,i,i);
    fprintf(outfile,"\n");
  }

  fprintf(outfile,"\n");

  // kernal of loop struct
  i=0;
  for(list<XMLString>::const_iterator pXMLString = mySamplingMomentsList.begin();pXMLString != mySamplingMomentsList.end(); pXMLString++) {
    for(j=0;j<myNDims;j++) {
	fprintf(outfile,"	");
    }
    fprintf(outfile,"%s=_%s_raw[_%s_raw_index_pointer+%li];\n",pXMLString->c_str(),mgFieldName,mgFieldName,i);
    i++;
  }
  fprintf(outfile,"\n");

  char indent[64];
  for(i=0;i<myNDims;i++) {
      indent[i]=0x09;
  }
  indent[myNDims]=0;

  const char* dotre="";
  if(complexRawVector) {
    dotre=".re";
  }

  i=0;
  if(myPostCode.length()==0) {
    for(list<XMLString>::const_iterator pXMLString = mySamplingMomentsList.begin(); pXMLString != mySamplingMomentsList.end(); pXMLString++) {
      if(simulation()->parameters()->nPaths>1) {
	fprintf(outfile,"%s_to_field[_%s_processed_index_pointer + %li] += %s%s;\n",
		indent,mgFieldName,i,pXMLString->c_str(),dotre);
	fprintf(outfile,"%s_to_field_sd[_%s_processed_index_pointer + %li] += %s%s*%s%s;\n",
		indent,mgFieldName,i,pXMLString->c_str(),dotre,pXMLString->c_str(),dotre);
      }
      else {
	fprintf(outfile,"%s_to_field[_%s_processed_index_pointer + %li] = %s%s;\n",
		indent,mgFieldName,i,pXMLString->c_str(),dotre);
      }
      i++;
    }
  }
  else {
    fprintf(outfile,"// ********** Post propagation code *************\n");
    fprintf(outfile,"%s\n",myPostCode.c_str());
    fprintf(outfile,"// **********************************************\n");
    for(list<XMLString>::const_iterator pXMLString = myPostMomentsList.begin(); pXMLString != myPostMomentsList.end(); pXMLString++) {
      if(simulation()->parameters()->nPaths>1) {
	fprintf(outfile,"%s_to_field[_%s_processed_index_pointer + %li] += %s;\n",indent,mgFieldName,i,pXMLString->c_str());
	fprintf(outfile,"%s_to_field_sd[_%s_processed_index_pointer + %li] += %s*%s;\n",indent,mgFieldName,i,pXMLString->c_str(),pXMLString->c_str());
      }
      else
	fprintf(outfile,"%s_to_field[_%s_processed_index_pointer + %li] = %s;\n",indent,mgFieldName,i,pXMLString-