xmdsfield.cc

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

/*
 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.31 2004/08/02 08:47:23 paultcochrane 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());
    }
  }
};

// ******************************************************************************