xmdsfield.cc

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

void xmdsField::outputSampleCount() const {
  if(debugFlag) {
    printf("xmdsField::outputSampleCount\n");
  }

  if(mySamplesList.size() != simulation()->output()->nMomentGroups()) {
    throw xmdsException("number of samples in field element does not match number of output moment groups");
  }

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

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

  return myNeedsFFTWPlans;
};

// ******************************************************************************
void xmdsField::writeInitialisationCalls(
					 FILE *const outfile,
					 const char *const indent) const {
  if(debugFlag) {
    printf("xmdsField::writeInitialisationCalls\n");
  }

  for(list<const xmdsVector*>::const_iterator ppxmdsVector = myVectorsList.begin(); ppxmdsVector != myVectorsList.end(); ppxmdsVector++) {
    (*ppxmdsVector)->writeInitialisationCall(outfile,indent);
  }

  fprintf(outfile,"\n");
};

// ******************************************************************************
void xmdsField::writeSampleCalls(
				 FILE *const outfile,
				 const char *const indent) const {
  if(debugFlag) {
    printf("xmdsField::writeSampleCalls\n");
  }

  unsigned long i=0;
  for(list<unsigned long>::const_iterator pULong = mySamplesList.begin(); pULong != mySamplesList.end(); pULong++) {
    if(*pULong>0) {
      fprintf(outfile,"%s_mg%li_sample();\n",indent,i);
    }
    i++;
  }

  fprintf(outfile,"\n");
};

// ******************************************************************************
void xmdsField::assignActiveVectorPointers(
					   FILE *const outfile,
					   const char *const tempVectorName) const {
  if(debugFlag) {
    printf("xmdsField::assignActiveVectorPointers\n");
  }

  for(list<const xmdsVector*>::const_iterator ppxmdsVector = myVectorsList.begin(); ppxmdsVector != myVectorsList.end(); ppxmdsVector++) {
    fprintf(outfile,"_active_%s_%s = %s_%s_%s;\n",
	    myName.c_str(),(*ppxmdsVector)->name()->c_str(),tempVectorName,myName.c_str(),(*ppxmdsVector)->name()->c_str());
  }

  fprintf(outfile,"\n");
};

// ******************************************************************************
const XMLString* xmdsField::name() const {
  if(debugFlag) {
    printf("xmdsField::name\n");
  }
	
  return &myName;
};

// ******************************************************************************
const xmdsFieldGeometry* xmdsField::geometry() const{
  if(debugFlag) {
    printf("xmdsField::geometry\n");
  }

  return &myGeometry;
};

// ******************************************************************************
bool xmdsField::getVector(
			  const XMLString& vectorName,
			  const xmdsVector*& theVector) const {
  if(debugFlag) {
    printf("xmdsField::getVector\n");
  }

  theVector=0;

  for(list<const xmdsVector*>::const_iterator ppxmdsVector = myVectorsList.begin(); ppxmdsVector != myVectorsList.end(); ppxmdsVector++) {
    if(*(*ppxmdsVector)->name()==vectorName) {
      theVector = *ppxmdsVector;
      return 1;
    }
  }
  return 0;
};

// ******************************************************************************
void xmdsField::vectorNames(
			    list<XMLString>& vectorNamesList) const {
  if(debugFlag) {
    printf("xmdsField::vectorNames\n");
  }

  vectorNamesList.clear();

  for(list<const xmdsVector*>::const_iterator ppxmdsVector = myVectorsList.begin(); ppxmdsVector != myVectorsList.end(); ppxmdsVector++) {
    vectorNamesList.push_back(*(*ppxmdsVector)->name());
  }
};

// ******************************************************************************
void xmdsField::vectors2space(
			      FILE *const outfile,
			      const unsigned long& space,
			      const list<XMLString>& vectorNamesList,
			      const char *const indent) const {
  if(debugFlag) {
    printf("xmdsField::vectors2space\n");
  }

  if(myGeometry.nDims()==0) {
    return;
  }

  for(list<XMLString>::const_iterator pXMLString = vectorNamesList.begin(); pXMLString != vectorNamesList.end(); pXMLString++) {

    const xmdsVector* nextVector;

    if(!getVector(*pXMLString,nextVector)) {
      throw xmdsException("Internal error in xmdsField::vectors2space: cannot find vector");
    }

    if(nextVector->needsFFTWRoutines()) {
      fprintf(outfile,"%s_%s_%s_go_space(%li);\n",indent,myName.c_str(),pXMLString->c_str(),space);
      fprintf(outfile,"\n");
    }
  }
};

// ******************************************************************************
void xmdsField::openLoops(
			  FILE *const outfile,
			  unsigned long space,
						  const list<XMLString>& vectorNamesList) const {
	if(debugFlag) {
		printf("xmdsField::openLoops\n");
	}
	
	for(list<XMLString>::const_iterator pXMLString = vectorNamesList.begin(); pXMLString != vectorNamesList.end(); pXMLString++) {
		fprintf(outfile,"unsigned long _%s_%s_index_pointer=0;\n",myName.c_str(),pXMLString->c_str());
	}
	fprintf(outfile,"\n");
	
	bool swapped=false;
	if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
		if((space&1)&((space>>1)&1)) {
			swapped=true;
		}
	}
	
	if((!simulation()->parameters()->stochastic)&(simulation()->parameters()->usempi)&swapped) {
		fprintf(outfile,
				"double k%s = local_y_start_after_transpose*_%s_dk1;\n"
				"if(k%s>((_%s_lattice1-1)/2 + 0.1)*_%s_dk1)\n"
				"	 k%s -= _%s_lattice1*_%s_dk1;\n"
				"\n"    
				"for(long _i0=0; _i0<local_ny_after_transpose; _i0++) {\n\n",
				myGeometry.dimension(1)->name.c_str(),myName.c_str(),
				myGeometry.dimension(1)->name.c_str(),myName.c_str(),myName.c_str(),
				myGeometry.dimension(1)->name.c_str(),myName.c_str(),myName.c_str());
		
		space >>= 1;
		
		fprintf(outfile,
				"	"
				"double k%s = 0;\n"
				"\n"
				"	"
				"for(long _i1=0; _i1<_%s_lattice0; _i1++) {\n"
				"\n",
				myGeometry.dimension(0)->name.c_str(),
				myName.c_str());
		
		space >>= 1;
		
		for(unsigned long i=2; i<myGeometry.nDims(); i++) {
			
			for(unsigned long j=0; j<i; j++) {
				fprintf(outfile,"	");
			}
			if(space&1) {
				fprintf(outfile,"double k%s = 0;\n",myGeometry.dimension(i)->name.c_str());
			}
			else {
				fprintf(outfile,"double %s = _%s_xmin%li;\n",myGeometry.dimension(i)->name.c_str(),myName.c_str(),i);
			}
			fprintf(outfile,"\n");
			
			for(unsigned long j=0; j<i; j++) {
				fprintf(outfile,"	");
			}
			fprintf(outfile,"for(long _i%li=0; _i%li<_%s_lattice%li; _i%li++) {\n",i,i,myName.c_str(),i,i);
			fprintf(outfile,"\n");
			
			space >>= 1;
		}
    }
	else if((!simulation()->parameters()->stochastic)&(simulation()->parameters()->usempi)) {
		if(space&1) {
			fprintf(outfile,"double k%s = local_x_start*_%s_dk0;\n",myGeometry.dimension(0)->name.c_str(),myName.c_str());
			fprintf(outfile,"if(k%s>((_%s_lattice0-1)/2 + 0.1)*_%s_dk0)\n",
					myGeometry.dimension(0)->name.c_str(),myName.c_str(),myName.c_str());	
			fprintf(outfile,"	k%s -= _%s_lattice0*_%s_dk0;\n",
					myGeometry.dimension(0)->name.c_str(),myName.c_str(),myName.c_str());
		}
		else
			fprintf(outfile,"double %s = _%s_xmin0+local_x_start*_%s_dx0;\n",myGeometry.dimension(0)->name.c_str(),myName.c_str(),myName.c_str());
		fprintf(outfile,"\n");
		
		fprintf(outfile,"for(long _i0=0; _i0<local_nx; _i0++) {\n");
		fprintf(outfile,"\n");
		
		space >>= 1;
		
		for(unsigned long i=1; i<myGeometry.nDims(); i++) {
			
			for(unsigned long j=0; j<i; j++) {
				fprintf(outfile,"	");
			}
			if(space&1) {
				fprintf(outfile,"double k%s = 0;\n",myGeometry.dimension(i)->name.c_str());
			}
			else {
				fprintf(outfile,"double %s = _%s_xmin%li;\n",myGeometry.dimension(i)->name.c_str(),myName.c_str(),i);
			}
			fprintf(outfile,"\n");
			
			for(unsigned long j=0; j<i; j++) {
				fprintf(outfile,"	");
			}
			fprintf(outfile,"for(long _i%li=0; _i%li<_%s_lattice%li; _i%li++) {\n",i,i,myName.c_str(),i,i);
			fprintf(outfile,"\n");
			
			space >>= 1;
		}
    }
	else {
		for(unsigned long i=0; i<myGeometry.nDims(); i++) {
			
			for(unsigned long j=0; j<i; j++) {
				fprintf(outfile,"	");
			}
			if(space&1) {
				fprintf(outfile,"double k%s = 0;\n",myGeometry.dimension(i)->name.c_str());
			}
			else {
				fprintf(outfile,"double %s = _%s_xmin%li;\n",myGeometry.dimension(i)->name.c_str(),myName.c_str(),i);
			}
			fprintf(outfile,"\n");
			
			for(unsigned long j=0; j<i; j++) {
				fprintf(outfile,"	");
			}
			fprintf(outfile,"for(long _i%li=0; _i%li<_%s_lattice%li; _i%li++) {\n",i,i,myName.c_str(),i,i);
			fprintf(outfile,"\n");
			
			space >>= 1;
		}
	}
						  };

// ******************************************************************************
void xmdsField::closeLoops(
			   FILE *const outfile,
			   unsigned long space,
			   const list<XMLString>& vectorNamesList) const {
  if(debugFlag) {
    printf("xmdsField::closeLoops\n");
  }

  if(myGeometry.nDims()>0) {
    fprintf(outfile,"\n");
    for(list<XMLString>::const_iterator pXMLString = vectorNamesList.begin(); pXMLString != vectorNamesList.end(); pXMLString++) {
      for(unsigned long j=0; j<myGeometry.nDims(); j++) {
	  fprintf(outfile,"	");
      }
      fprintf(outfile,"_%s_%s_index_pointer += _%s_%s_ncomponents;\n",
	      myName.c_str(),pXMLString->c_str(),myName.c_str(),pXMLString->c_str());
    }
  }

  bool *isKSpace = new bool[myGeometry.nDims()];
  unsigned long tempSpace = space;
  for(long unsigned int i=0; i<myGeometry.nDims(); i++) {
	 if(tempSpace&1) {
		isKSpace[i] = true;
		}
	 else {
		isKSpace[i] = false;
		}
	tempSpace>>=1;
  }

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

  if((!simulation()->parameters()->stochastic)&(simulation()->parameters()->usempi)&swapped) {
      for(unsigned long i=myGeometry.nDims(); i>2; i--) {
		  fprintf(outfile,"\n");
		  if(isKSpace[i-1]) {
			  for(unsigned long j=0; j<i; j++) {
				  fprintf(outfile,"	");
			  }
			  fprintf(outfile,"k%s += _%s_dk%li;\n",myGeometry.dimension(i-1)->name.c_str(),myName.c_str(),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",
					  myGeometry.dimension(i-1)->name.c_str(),myName.c_str(),i-1,myName.c_str(),i-1);	
			  
			  for(unsigned long j=0; j<i; j++) {
				  fprintf(outfile,"	");
			  }
			  fprintf(outfile,"	k%s -= _%s_lattice%li*_%s_dk%li;\n",
					  myGeometry.dimension(i-1)->name.c_str(),myName.c_str(),i-1,myName.c_str(),i-1);
		  }
		  else {
			  for(unsigned long j=0; j<i; j++) {
				  fprintf(outfile,"	");
			  }
			  fprintf(outfile,"%s += _%s_dx%li;\n",myGeometry.dimension(i-1)->name.c_str(),myName.c_str(),i-1);
		  }
		  
		  for(unsigned long j=0; j<i; j++) {
			  fprintf(outfile,"	");
		  }
		  fprintf(outfile,"}\n");
	  }
	  
      fprintf(outfile,"\n");
	  
      for(unsigned long j=0; j<2; j++) {
		  fprintf(outfile,"	");
      }
      fprintf(outfile,"k%s += _%s_dk0;\n",myGeometry.dimension(0)->name.c_str(),myName.c_str());
	  
      for(unsigned long j=0; j<2; j++) {
		  fprintf(outfile,"	");
      }
      fprintf(outfile,"if(k%s>((_%s_lattice0-1)/2 + 0.1)*_%s_dk0)\n",
			  myGeometry.dimension(0)->name.c_str(),myName.c_str(),myName.c_str());	
	  
      for(unsigned long j=0; j<2; j++) {
		  fprintf(outfile,"	");
      }
      fprintf(outfile,"	k%s -= _%s_lattice0*_%s_dk0;\n",
			  myGeometry.dimension(0)->name.c_str(),myName.c_str(),myName.c_str());
	  
      for(unsigned long j=0; j<2; j++) {
		  fprintf(outfile,"	");
      }
      fprintf(outfile,"}\n");
	  
      fprintf(outfile,"\n");
	  
      fprintf(outfile,"	");
      fprintf(outfile,"k%s += _%s_dk1;\n",myGeometry.dimension(1)->name.c_str(),myName.c_str());
	  
      fprintf(outfile,"	");
      fprintf(outfile,"if(k%s>((_%s_lattice1-1)/2 + 0.1)*_%s_dk1)\n",
			  myGeometry.dimension(1)->name.c_str(),myName.c_str(),myName.c_str());	
	  
      fprintf(outfile,"	");
      fprintf(outfile,"	k%s -= _%s_lattice1*_%s_dk1;\n",
			  myGeometry.dimension(1)->name.c_str(),myName.c_str(),myName.c_str());
	  
      fprintf(outfile,"	");
      fprintf(outfile,"}\n");
  }
  else
    {
      for(unsigned long i=myGeometry.nDims(); i>0; i--) {
		  
		  fprintf(outfile,"\n");