/*
Copyright (C) 2000-2004
Code contributed by Greg Collecutt, Joseph Hope and Paul Cochrane
This file is part of xmds.
This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 2
of the License, or (at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
*/
/*
$Id: xmdsfield.cc,v 1.32 2005/04/27 05:44:19 joehope Exp $
*/
/*! @file xmdsfield.cc
@brief Field element parsing classes and methods
More detailed explanation...
*/
#include<xmlbasics.h>
#include<dom3.h>
#include<xmdsutils.h>
#include<xmdsclasses.h>
// ******************************************************************************
// ******************************************************************************
// xmdsFieldGeometry
// ******************************************************************************
// ******************************************************************************
extern bool debugFlag;
// ******************************************************************************
void xmdsFieldGeometry::clear() {
if(debugFlag) {
printf("xmdsFieldGeometry::clear\n");
}
myDimensionsList.clear();
};
// ******************************************************************************
unsigned long xmdsFieldGeometry::nDims() const {
if(debugFlag) {
printf("xmdsFieldGeometry::nDims\n");
}
return myDimensionsList.size();
};
// ******************************************************************************
const dimensionStruct* xmdsFieldGeometry::dimension(
const unsigned long& index) const {
if(debugFlag) {
printf("xmdsFieldGeometry::dimension\n");
}
if(index>=myDimensionsList.size()) {
throw xmdsException("Internal range error in xmdsFieldGeometry::dimension");
}
list<dimensionStruct>::const_iterator pdimensionStruct = myDimensionsList.begin();
for(unsigned long i=0; i<index; i++) {
pdimensionStruct++;
}
return &*pdimensionStruct;
};
// ******************************************************************************
void xmdsFieldGeometry::setDimension(
const unsigned long& index,
const dimensionStruct& newDimension) {
if(debugFlag) {
printf("xmdsFieldGeometry::setDimension\n");
}
if(index>myDimensionsList.size()) {
throw xmdsException("Internal range error in xmdsFieldGeometry::setDimension");
}
if(index==myDimensionsList.size()) {
myDimensionsList.push_back(newDimension);
}
list<dimensionStruct>::iterator pdimensionStruct = myDimensionsList.begin();
for(unsigned long i=0; i<index; i++) {
pdimensionStruct++;
}
*pdimensionStruct=newDimension;
};
// ******************************************************************************
void xmdsFieldGeometry::addDimension(
const dimensionStruct& newDimension) {
if(debugFlag) {
printf("xmdsFieldGeometry::addDimension\n");
}
myDimensionsList.push_back(newDimension);
};
// ******************************************************************************
bool xmdsFieldGeometry::getDimNumber(
const XMLString& dimName,
unsigned long& dimNumber) const {
if(debugFlag) {
printf("xmdsFieldGeometry::getDimNumber\n");
}
unsigned long i = 0;
for(list<dimensionStruct>::const_iterator pDim = myDimensionsList.begin(); pDim != myDimensionsList.end(); pDim++) {
if(pDim->name == dimName) {
dimNumber=i;
return 1;
}
i++;
}
return 0;
};
// ******************************************************************************
unsigned long xmdsFieldGeometry::fullSpace() const {
if(debugFlag) {
printf("xmdsFieldGeometry::fullSpace\n");
}
unsigned long two2n=1;
for(unsigned long i=0; i<myDimensionsList.size(); i++) {
two2n *= 2;
}
return two2n-1;
};
// ******************************************************************************
// ******************************************************************************
// xmdsField Public
// ******************************************************************************
// ******************************************************************************
long nxmdsFields=0; //!< Number of xmds field objects
// ******************************************************************************
xmdsField::xmdsField(
const xmdsSimulation *const yourSimulation,
const bool& yourVerboseMode) :
xmdsElement(yourSimulation,yourVerboseMode) {
if(debugFlag) {
nxmdsFields++;
printf("xmdsField::xmdsField\n");
printf("nxmdsFields=%li\n",nxmdsFields);
}
};
// ******************************************************************************
xmdsField::~xmdsField() {
if(debugFlag) {
nxmdsFields--;
printf("xmdsField::~xmdsField\n");
printf("nxmdsFields=%li\n",nxmdsFields);
}
// destroy xmdsVectors, since they are not managed by the xmdsElement class
for(list<const xmdsVector*>::const_iterator ppxmdsVector = myVectorsList.begin(); ppxmdsVector != myVectorsList.end(); ppxmdsVector++) {
delete (*ppxmdsVector);
}
};
// ******************************************************************************
void xmdsField::processElement(
const Element *const yourElement) {
if(debugFlag) {
printf("xmdsField::processElement\n");
}
if(verbose()) {
printf("Processing field element ...\n");
}
myNeedsFFTWPlans=0;
list<XMLString> anXMLStringList;
list<XMLString>::const_iterator pXMLString;
// ************************************
// find name
getAssignmentStrings(yourElement,"name",0,1,anXMLStringList);
if(anXMLStringList.size()==1) {
myName = *anXMLStringList.begin();
if(verbose()) {
printf("field name is '%s'\n",myName.c_str());
}
}
else {
myName = "main";
printf("field name defaulting to '%s'\n",myName.c_str());
}
list<unsigned long> aLatticeList;
list<unsigned long>::const_iterator pULong;
list<domainStruct> aDomainList;
list<domainStruct>::const_iterator pDomain;
// ************************************
// find transverse dimensions
getAssignmentStrings(yourElement,"dimensions",0,0,anXMLStringList);
if(anXMLStringList.size()>63) {
throw xmdsException(yourElement,"xmds currently limited to 63 transverse dimensions");
}
if(anXMLStringList.size()>0) {
unsigned long i;
// ************************************
// find transverse lattice
getAssignmentULongs(yourElement,"lattice",1,anXMLStringList.size(),aLatticeList);
i=0;
for(pULong = aLatticeList.begin(); pULong != aLatticeList.end(); pULong++) {
if(*pULong < 2) {
sprintf(errorMessage(), "lattice for dimension #%li must be > 1",i+1);
throw xmdsException(yourElement,errorMessage());
}
if(verbose()) {
printf("transverse lattice #%li has %li points\n",i+1,*pULong);
}
i++;
}
// ************************************
// find transverse domains
getAssignmentDomains(yourElement,"domains",1,anXMLStringList.size(),aDomainList);
if(verbose()) {
i=0;
for(pDomain = aDomainList.begin(); pDomain != aDomainList.end(); pDomain++) {
printf("transverse domain #%li:\n",i+1);
printf(" begin = %f\n",pDomain->begin);
printf(" end = %f\n",pDomain->end);
i++;
}
}
}
else {
if(verbose()) {
printf("No transverse dimensions found.\n");
}
}
// build field geometry dimensions
myGeometry.clear();
pXMLString = anXMLStringList.begin();
pULong = aLatticeList.begin();
pDomain = aDomainList.begin();
for(unsigned long i=0; i< anXMLStringList.size(); i++) {
dimensionStruct nextDimension;
for(unsigned long j=0;j<i;j++) {
if(*pXMLString==myGeometry.dimension(j)->name) {
sprintf(errorMessage(),"duplicate dimension '%s'",pXMLString->c_str());
throw xmdsException(yourElement,errorMessage());
}
}
if(verbose()) {
printf("adding transverse dimension '%s'\n",pXMLString->c_str());
}
nextDimension.name = *pXMLString;
nextDimension.lattice = *pULong;
nextDimension.domain = *pDomain;
myGeometry.addDimension(nextDimension);
pXMLString++;
pULong++;
pDomain++;
}
// ************************************
// find samples
getAssignmentULongs(yourElement,"samples",1,0,mySamplesList);
// ************************************
// find and process vector elements
// ************************************
const NodeList* candidateElements;
candidateElements = yourElement->getElementsByTagName("vector",0);
if(candidateElements->length()==0) {
throw xmdsException(yourElement,"at least one vector expected");
}
for(unsigned long i=0; i<candidateElements->length(); i++) {
const Element* nextElement = dynamic_cast<const Element*>(candidateElements->item(i));
xmdsElement* newxmdsVectorElement = createxmdsVectorElement();
newxmdsVectorElement->processElement(nextElement);
}
};
// ******************************************************************************
void xmdsField::processVectors(
const list<XMLString>& vectorNamesList,
const unsigned long& space) const {
if(debugFlag) {
printf("xmdsField::processVectors\n");
}
for(list<XMLString>::const_iterator pXMLString = vectorNamesList.begin(); pXMLString != vectorNamesList.end(); pXMLString++) {
const xmdsVector* nextVector;
if(!getVector(*pXMLString,nextVector)) {
sprintf(errorMessage(),"Vector '%s' unknown",pXMLString->c_str());
throw xmdsException(errorMessage());
}
for(list<XMLString>::const_iterator pXMLString2 = vectorNamesList.begin(); pXMLString2 != pXMLString; pXMLString2++) {
if(*pXMLString == *pXMLString2) {
sprintf(errorMessage(),"Duplicate vector '%s'",pXMLString->c_str());
throw xmdsException(errorMessage());
}
}
if(nextVector->initialSpace() != space){
if(nextVector->vectorType() != COMPLEX) {
const char* spaceDescription;
if(space == 0) {
spaceDescription = "x-space";
}
else if(space == myGeometry.fullSpace()) {
spaceDescription = "k-space";
}
else {
spaceDescription = "mixed fourier space";
}
sprintf(errorMessage(),"Cannot convert vector '%s' to %s since it is not of type 'complex'",
pXMLString->c_str(),spaceDescription);
throw xmdsException(errorMessage());
}
myNeedsFFTWPlans = 1;
nextVector->setNeedsFFTWRoutines();
}
// Always write FFTW Plans for main vector
// (otherwise a rare exception for deterministic MPI which uses FFTW to split field across processes)
// Could make this MPI specifc, or else define the appropriate
if(pXMLString->c_str()=="main")
myNeedsFFTWPlans = 1;
if(verbose()) {
printf("adding vector '%s' to accessible vectors list\n",pXMLString->c_str());
}
}
};
// ******************************************************************************
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");
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");
}
}
delete isKSpace;
};
// ******************************************************************************
void xmdsField::writePlanCreationCalls(
FILE *const outfile,
const bool& useFFTWMeasure,
const bool& useWisdom) const {
if(debugFlag) {
printf("xmdsField::writePlanCreationCalls\n");
}
if(!myNeedsFFTWPlans) {
return;
}
const char *useMeasureStr = "";
const char *useWisdomStr = "";
//const char *wisdomFname = simulation()->parameters()->simulationName.c_str();
if(useFFTWMeasure) {
useMeasureStr = "|FFTW_MEASURE";
}
if(useWisdom) {
useWisdomStr = "|FFTW_USE_WISDOM";
}
for(unsigned long i=0;i<myGeometry.nDims();i++) {
fprintf(outfile,"_fftw_lattice[%li]=_%s_lattice%li;\n",i,myName.c_str(),i);
}
if((!(simulation()->parameters()->stochastic))&(simulation()->parameters()->usempi)){
if (useWisdom) {
fprintf(outfile,
"{\n"
"\nFILE *fWisdom;\n"
"ifstream fIn;\n"
"ofstream fOut;\n"
"// this is a dodgy hack, but it should work\n"
"string findHostString, findHomeString, hostStuff, homeStuff, rmString, rankString;\n"
"ostringstream outString;\n"
"outString << rank;\n"
"rankString = outString.str();\n"
"hostStuff = \"host_rank_\" + rankString + \".stuff\";\n"
"findHostString = \"uname -n > \" + hostStuff;\n"
"system(findHostString.c_str());\n"
"homeStuff = \"home_rank_\" + rankString + \".stuff\";\n"
"findHomeString = \"echo $HOME > \" + homeStuff;\n"
"system(findHomeString.c_str());\n"
"string hostName, homeDir;\n"
"fIn.open(hostStuff.c_str());\n"
"if (fIn.fail()) {\n"
" // do something\n"
"}\n"
"fIn >> hostName;\n"
"fIn.close();\n"
"fIn.open(homeStuff.c_str());\n"
"if (fIn.fail()) {\n"
" // do something\n"
"}\n"
"fIn >> homeDir;\n"
"fIn.close();\n"
"rmString = \"rm \" + hostStuff + \" \" + homeStuff;\n"
"system(rmString.c_str());\n\n"
"string testFname, wisdomFname;\n"
"testFname = homeDir + \"/.xmds/wisdom/test-\" + rankString;\n"
"fOut.open(testFname.c_str());\n"
"if (fOut.fail()) {\n"
" cout << \"Warning: ~/.xmds/wisdom directory doesn't seem to exist.\\n\";\n"
" cout << \"Using current directory instead\\n\";\n"
" wisdomFname = hostName + \".wisdom\";\n"
"}\n"
"else {\n"
" fOut.close();\n"
" rmString = \"rm \" + testFname;\n"
" system(rmString.c_str());\n"
" wisdomFname = homeDir + \"/.xmds/wisdom/\" + hostName + \".wisdom\";\n"
"}\n"
"printf(\"Performing fftw calculations\\n\");\n"
"if ((fWisdom = fopen(wisdomFname.c_str(), \"r\")) != NULL) {\n"
" cout << \"Standing upon the shoulders of giants... (Importing wisdom)\\n\";\n"
" fftw_import_wisdom_from_file(fWisdom);\n"
" fclose(fWisdom);\n"
"}\n"
);
}
fprintf(outfile,"printf(\"Making forward plan\\n\");\n");
fprintf(outfile,"_%s_forward_plan = fftwnd_mpi_create_plan(MPI_COMM_WORLD,_%s_ndims,_fftw_lattice,FFTW_FORWARD,FFTW_IN_PLACE%s%s);\n",myName.c_str(),myName.c_str(),useMeasureStr,useWisdomStr);
// Using FFTW_TRANSPOSED_ORDER, so backward plan has different array sizes
fprintf(outfile,"_fftw_lattice[0]=_%s_lattice1;\n",myName.c_str());
fprintf(outfile,"_fftw_lattice[1]=_%s_lattice0;\n",myName.c_str());
for(unsigned long i=2;i<myGeometry.nDims();i++) {
fprintf(outfile,"_fftw_lattice[%li]=_%s_lattice%li;\n",i,myName.c_str(),i);
}
fprintf(outfile,"printf(\"Making backward plan\\n\");\n");
fprintf(outfile,"_%s_backward_plan = fftwnd_mpi_create_plan(MPI_COMM_WORLD,_%s_ndims,_fftw_lattice,FFTW_BACKWARD,FFTW_IN_PLACE%s%s);\n",myName.c_str(),myName.c_str(),useMeasureStr,useWisdomStr);
if (useWisdom) {
fprintf(outfile,"printf(\"Keeping accumulated wisdom\\n\");\n");
fprintf(outfile,"fWisdom = fopen(wisdomFname.c_str(), \"w\");\n");
fprintf(outfile,"fftw_export_wisdom_to_file(fWisdom);\n");
fprintf(outfile,"fclose(fWisdom);\n");
fprintf(outfile,"}\n");
fprintf(outfile,"printf(\"Finished fftw calculations\\n\");\n");
}
}
else {
if (useWisdom) {
fprintf(outfile,
"{\n"
"FILE *fWisdom;\n"
"ifstream fIn;\n"
"ofstream fOut;\n"
"// this is a dodgy hack, but it should work\n"
"string findHostString, findHomeString, hostStuff, homeStuff, rmString, rankString;\n");
if (!simulation()->parameters()->usempi) {
fprintf(outfile, "int rank = 0;\n");
}
fprintf(outfile,
"ostringstream outString;\n"
"outString << rank;\n"
"rankString = outString.str();\n"
"hostStuff = \"host_rank_\" + rankString + \".stuff\";\n"
"findHostString = \"uname -n > \" + hostStuff;\n"
"system(findHostString.c_str());\n"
"homeStuff = \"home_rank_\" + rankString + \".stuff\";\n"
"findHomeString = \"echo $HOME > \" + homeStuff;\n"
"system(findHomeString.c_str());\n"
"string hostName, homeDir;\n"
"fIn.open(hostStuff.c_str());\n"
"if (fIn.fail()) {\n"
" // do something\n"
"}\n"
"fIn >> hostName;\n"
"fIn.close();\n"
"fIn.open(homeStuff.c_str());\n"
"if (fIn.fail()) {\n"
" // do something\n"
"}\n"
"fIn >> homeDir;\n"
"fIn.close();\n"
"rmString = \"rm \" + hostStuff + \" \" + homeStuff;\n"
"system(rmString.c_str());\n\n"
"string testFname, wisdomFname;\n"
"testFname = homeDir + \"/.xmds/wisdom/test-\" + rankString;\n"
"fOut.open(testFname.c_str());\n"
"if (fOut.fail()) {\n"
" cout << \"Warning: ~/.xmds/wisdom directory doesn't seem to exist.\\n\";\n"
" cout << \"Using current directory instead\\n\";\n"
" wisdomFname = hostName + \".wisdom\";\n"
"}\n"
"else {\n"
" fOut.close();\n"
" wisdomFname = homeDir + \"/.xmds/wisdom/\" + hostName + \".wisdom\";\n"
" rmString = \"rm \" + testFname;\n"
" system(rmString.c_str());\n"
"}\n"
"printf(\"Performing fftw calculations\\n\");\n"
"if ( (fWisdom = fopen(wisdomFname.c_str(), \"r\")) != NULL) {\n"
" cout << \"Standing upon the shoulders of giants... (Importing wisdom)\\n\";\n"
" fftw_import_wisdom_from_file(fWisdom);\n"
" fclose(fWisdom);\n"
"}\n");
}
fprintf(outfile,"printf(\"Making forward plan\\n\");\n");
fprintf(outfile,"_%s_forward_plan = fftwnd_create_plan(_%s_ndims,_fftw_lattice,FFTW_FORWARD,FFTW_IN_PLACE%s%s);\n",myName.c_str(),myName.c_str(),useMeasureStr,useWisdomStr);
fprintf(outfile,"printf(\"Making backward plan\\n\");\n");
fprintf(outfile,"_%s_backward_plan = fftwnd_create_plan(_%s_ndims,_fftw_lattice,FFTW_BACKWARD,FFTW_IN_PLACE%s%s);\n",myName.c_str(),myName.c_str(),useMeasureStr,useWisdomStr);
if (useWisdom) {
fprintf(outfile,
"printf(\"Keeping accumulated wisdom\\n\");\n"
"fWisdom = fopen(wisdomFname.c_str(),\"w\");\n"
"fftw_export_wisdom_to_file(fWisdom);\n"
"fclose(fWisdom);\n"
"}\n"
"printf(\"Finished fftw calculations\\n\");\n");
}
}
fprintf(outfile,"\n");
};
// ******************************************************************************
void xmdsField::writePlanDeletionCalls(
FILE *const outfile) const {
if(debugFlag) {
printf("xmdsField::writePlanDeletionCalls\n");
}
if(!myNeedsFFTWPlans) {
return;
}
if((!(simulation()->parameters()->stochastic))&(simulation()->parameters()->usempi)){
fprintf(outfile,"fftwnd_mpi_destroy_plan(_%s_forward_plan);\n",myName.c_str());
fprintf(outfile,"fftwnd_mpi_destroy_plan(_%s_backward_plan);\n",myName.c_str());
}
else {
fprintf(outfile,"fftwnd_destroy_plan(_%s_forward_plan);\n",myName.c_str());
fprintf(outfile,"fftwnd_destroy_plan(_%s_backward_plan);\n",myName.c_str());
}
fprintf(outfile,"\n");
};
// ******************************************************************************
xmdsVector* xmdsField::createxmdsVector() {
if(debugFlag) {
printf("xmdsField::createxmdsVector\n");
}
xmdsVector* newVector = new xmdsVector(this);
myVectorsList.push_back(newVector);
return newVector;
};
// ******************************************************************************
// ******************************************************************************
// xmdsField protected
// ******************************************************************************
// ******************************************************************************
// ******************************************************************************
void xmdsField::writeDefines(
FILE *const outfile) const {
if(debugFlag) {
printf("xmdsField::writeDefines\n");
}
fprintf(outfile,"// ********************************************************\n");
fprintf(outfile,"// field %s defines\n",myName.c_str());
fprintf(outfile,"\n");
fprintf(outfile,"#define _%s_ndims %li\n",myName.c_str(),myGeometry.nDims());
for(unsigned long i=0; i<myGeometry.nDims(); i++) {
const dimensionStruct* dimI = myGeometry.dimension(i);
fprintf(outfile,"#define _%s_lattice%li %li\n",myName.c_str(),i,dimI->lattice);
fprintf(outfile,"#define _%s_xmin%li %e\n",myName.c_str(),i,dimI->domain.begin);
fprintf(outfile,"#define _%s_dx%li ((%e - %e)/(double)_%s_lattice%li)\n",myName.c_str(),i,dimI->domain.end,dimI->domain.begin,myName.c_str(),i);
fprintf(outfile,"#define _%s_dk%li (2*M_PI/(%e - %e))\n",myName.c_str(),i,dimI->domain.end,dimI->domain.begin);
fprintf(outfile,"#define d%s _%s_dx%li\n",dimI->name.c_str(),myName.c_str(),i);
fprintf(outfile,"#define dk%s _%s_dk%li\n",dimI->name.c_str(),myName.c_str(),i);
}
fprintf(outfile,"\n");
for(list<const xmdsVector*>::const_iterator ppxmdsVector = myVectorsList.begin(); ppxmdsVector != myVectorsList.end(); ppxmdsVector++) {
(*ppxmdsVector)->writeDefines(outfile);
}
fprintf(outfile,"\n");
};
// ******************************************************************************
void xmdsField::writeGlobals(
FILE *const outfile) const {
if(debugFlag) {
printf("xmdsField::writeGlobals\n");
}
fprintf(outfile,"// ********************************************************\n");
fprintf(outfile,"// field %s globals\n",myName.c_str());
fprintf(outfile,"\n");
if(myGeometry.nDims()>0) {
fprintf(outfile,"const unsigned long _%s_size=_%s_lattice0",myName.c_str(),myName.c_str());
for(unsigned long i=1; i<myGeometry.nDims(); i++) {
fprintf(outfile,"*_%s_lattice%li",myName.c_str(),i);
}
fprintf(outfile,";\n");
}
else
fprintf(outfile,"const unsigned long _%s_size=1;\n",myName.c_str());
fprintf(outfile,"\n");
if(myNeedsFFTWPlans) {
if((!(simulation()->parameters()->stochastic))&(simulation()->parameters()->usempi)) {
fprintf(outfile,"fftwnd_mpi_plan _%s_forward_plan;\n",myName.c_str());
fprintf(outfile,"fftwnd_mpi_plan _%s_backward_plan;\n",myName.c_str());
}
else {
fprintf(outfile,"fftwnd_plan _%s_forward_plan;\n",myName.c_str());
fprintf(outfile,"fftwnd_plan _%s_backward_plan;\n",myName.c_str());
}
}
fprintf(outfile,"\n");
for(list<const xmdsVector*>::const_iterator ppxmdsVector = myVectorsList.begin(); ppxmdsVector != myVectorsList.end(); ppxmdsVector++) {
(*ppxmdsVector)->writeGlobals(outfile);
}
fprintf(outfile,"\n");
};
// ******************************************************************************
void xmdsField::writePrototypes(
FILE *const outfile) const {
if(debugFlag) {
printf("xmdsField::writePrototypes\n");
}
fprintf(outfile,"// ********************************************************\n");
fprintf(outfile,"// field %s prototypes\n",myName.c_str());
fprintf(outfile,"\n");
for(list<const xmdsVector*>::const_iterator ppxmdsVector = myVectorsList.begin(); ppxmdsVector != myVectorsList.end(); ppxmdsVector++) {
(*ppxmdsVector)->writePrototypes(outfile);
}
fprintf(outfile,"\n");
};
// ******************************************************************************
void xmdsField::writeRoutines(
FILE *const outfile) const {
if(debugFlag) {
printf("xmdsField::writeRoutines\n");
}
fprintf(outfile,"// ********************************************************\n");
fprintf(outfile,"// field %s routines\n",myName.c_str());
fprintf(outfile,"// ********************************************************\n");
fprintf(outfile,"\n");
for(list<const xmdsVector*>::const_iterator ppxmdsVector = myVectorsList.begin(); ppxmdsVector != myVectorsList.end(); ppxmdsVector++) {
(*ppxmdsVector)->writeRoutines(outfile);
}
fprintf(outfile,"\n");
};
// ******************************************************************************
void xmdsField::setName(
const XMLString& yourName) {
if(debugFlag) {
printf("xmdsField::setName\n");
}
myName=yourName;
};
// ******************************************************************************
void xmdsField::setGeometry(
const xmdsFieldGeometry& yourGeometry) {
if(debugFlag) {
printf("xmdsField::setGeometry\n");
}
myGeometry = yourGeometry;
};
// ******************************************************************************
// ******************************************************************************
// xmdsField private
// ******************************************************************************
// ******************************************************************************
// ******************************************************************************
xmdsVectorElement* xmdsField::createxmdsVectorElement() {
if( debugFlag) {
printf("xmdsField::createxmdsVectorElement\n");
}
xmdsVectorElement* newVectorElement = new xmdsVectorElement(simulation(),verbose(),this);
myVectorsList.push_back(newVectorElement);
return newVectorElement;
};
syntax highlighted by Code2HTML, v. 0.9.1