xmdsvectorelement.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: xmdsvectorelement.cc,v 1.17 2004/10/21 09:51:22 paultcochrane Exp $
*/

/*! @file xmdsvectorelement.cc
  @brief Vector element parsing classes and methods

  More detailed explanation...
*/

#include<xmlbasics.h>
#include<dom3.h>
#include<xmdsutils.h>
#include<xmdsclasses.h>

// ******************************************************************************
// ******************************************************************************
//                              xmdsVectorElement public
// ******************************************************************************
// ******************************************************************************

extern bool debugFlag;

long nxmdsVectorElements=0;   //!< The number of xmds vector elements

enum {
  ASCII = 1,
  BINARY = 2
};

// ******************************************************************************
xmdsVectorElement::xmdsVectorElement(
				     const xmdsSimulation *const yourSimulation,
				     const bool& yourVerboseMode,
				     const xmdsField *const yourField) :
  xmdsVector(yourField),
  xmdsElement(yourSimulation,yourVerboseMode) {
  if(debugFlag) {
    nxmdsVectorElements++;
    printf("xmdsVectorElement::xmdsVectorElement\n");
    printf("nxmdsVectorElements=%li\n",nxmdsVectorElements);
  }
};

// ******************************************************************************
xmdsVectorElement::~xmdsVectorElement() {
  if(debugFlag) {
    nxmdsVectorElements--;
    printf("xmdsVectorElement::~xmdsVectorElement\n");
    printf("nxmdsVectorElements=%li\n",nxmdsVectorElements);
  }
};

// ******************************************************************************
void xmdsVectorElement::processElement(
				       const Element *const yourElement) {
  if(debugFlag) {
    printf("xmdsVectorElement::processElement\n");
  }

  list<XMLString> anXMLStringList;
  XMLString anXMLString;

  if(verbose()) {
    printf("Processing vector element ...\n");
  }

  // ************************************
  // find name

  getAssignmentStrings(yourElement,"name",1,1,anXMLStringList);

  const xmdsVector* possibleTwin;

  if(field()->getVector(*anXMLStringList.begin(),possibleTwin)) {
    sprintf(errorMessage(),"a vector of name '%s' already exists in this field",anXMLStringList.begin()->c_str());
    throw xmdsException(yourElement,errorMessage());
  }

  setName(*anXMLStringList.begin());
  if(verbose()) {
    printf("vector name is '%s'\n",name()->c_str());
  }

  // ************************************
  // find type

  getAssignmentStrings(yourElement,"type",0,1,anXMLStringList);

  if(anXMLStringList.size()==1) {

    if(*anXMLStringList.begin()=="complex") {
      setVectorType(COMPLEX);
      if(verbose()) {
	printf("vector type is 'complex'\n");
      }
    }
    else if(*anXMLStringList.begin()=="double") {
      setVectorType(DOUBLE);
      if(verbose()) {
	printf("vector type is 'double'\n");
      }
    }
    else {
    }
  }
  else {
    setVectorType(COMPLEX);
    printf("vector type defaulting to 'complex'\n");
  }

  // ************************************
  // find components

  getAssignmentStrings(yourElement,"components",1,0,anXMLStringList);

  if(anXMLStringList.size()==0) {
    throw xmdsException(yourElement,"No vector components specified!");
  }

  setComponents(anXMLStringList);

  if(verbose()) {
    for(unsigned long i=0; i<nComponents(); i++) {
      printf("adding vector component '%s'\n",componentName(i)->c_str());
    }
  }

  // ************************************
  // find 'filename'

  getAssignmentStrings(yourElement,"filename",NOT_REQD,1,anXMLStringList);
  
  if(anXMLStringList.size()==1) {
    // we have a filename assignent

    getAttributeStrings(yourElement,"filename","format",NOT_REQD,anXMLString);    
    myInputFileFormat = ASCII;  // the default value

    if (anXMLString != EMPTY_STRING) {
    // a format attribute has been set
      if (anXMLString == "binary") {
	myInputFileFormat = BINARY;
      }
      else if (anXMLString == "bin") {
	myInputFileFormat = BINARY;
      }
      else if (anXMLString == "ascii") {
	myInputFileFormat = ASCII;
      }
      else if (anXMLString == "text") {
	myInputFileFormat = ASCII;
      }
      else if (anXMLString == "txt") {
	myInputFileFormat = ASCII;
      }
      else {
	throw xmdsException(yourElement, 
			    "Warning: Unknown file format attribute value in 'filename' tag\nI expected either 'binary' or 'ascii'.");
      }
    }
    
    myFileName=*anXMLStringList.begin();
    if(verbose()) {
      printf("vector initialisation from file '%s'\n",myFileName.c_str());
    }
    setInitialSpace(0);
  }
  else {
    // initialisation from code
    if(verbose()) {
      printf("initialisation from user code\n");
    }

    setInitialSpace(0);

    if(field()->geometry()->nDims()>0) {

      // ************************************
      // find space

      list<bool> aSpaceList;

      getAssignmentBools(yourElement,"fourier_space",0,field()->geometry()->nDims(),aSpaceList);

      if(aSpaceList.size() == 0) {
	printf("Initialisation space for vector '%s' defaulting to x-space.\n",name()->c_str());
      }
      else {
	list<bool>::const_iterator pBool = aSpaceList.begin();
	for(unsigned long i=0;i<field()->geometry()->nDims();i++) {
	  if(verbose()) {
	    if(*pBool) {
	      printf("initialisation will be performed with dimension #%li in fourier space\n",i+1);
	    }
	    else {
	      printf("initialisation will be performed with dimension #%li in normal space\n",i+1);
	    }
	  }
	  pBool++;
	}
	setInitialSpace(spaceList2ULong(aSpaceList));
      }
    }

    // ************************************
    // find vectors

    getAssignmentStrings(yourElement,"vectors",0,0,myVectorNamesList);

    field()->processVectors(myVectorNamesList,initialSpace());

    // ************************************
    // find code

    myCode=*yourElement->textContent(0);

    // check for non-white space charaters in code:

    if(myCode.isAllWhiteSpace()) {
      throw xmdsException(yourElement,"No initialisation code defined!");
    }

    if(verbose()) {
      printf("initialisation code loaded\n");
    }
  }
};

// ******************************************************************************
// ******************************************************************************
//                              xmdsVectorElement private
// ******************************************************************************
// ******************************************************************************

// ******************************************************************************
void xmdsVectorElement::writeInitialiseRoutine(
					       FILE *const outfile) const {
  if(debugFlag) {
    printf("xmdsVectorElement::writeInitialise\n");
  }

  const unsigned long nDims = field()->geometry()->nDims();
  const char *const fieldName = field()->name()->c_str();
  const char *const vectorName = name()->c_str();

  fprintf(outfile,
	  "// *************************\n"
	  "void _%s_%s_initialise() {\n\n",
	  fieldName,
	  vectorName);

  if(myFileName.length()>0) {
    // initialisation from file

    if(myInputFileFormat == ASCII) {

      fprintf(outfile,
	      "FILE* infile = fopen(\"%s\",\"r\");\n"
	      "\n"
	      "if (infile==0) {\n"
	      "\t printf(\"Error opening file '%s' for initialisation of vector '%s'\\n\");\n"
	      "\t exit(255);\n"
	      "\t }\n"
	      "\n"
	      "unsigned long _i0=0;\n"
	      "\n"
	      "while(_i0<_%s_size*_%s_%s_ncomponents) {\n",
	      myFileName.c_str(),
	      myFileName.c_str(),
	      vectorName,
	      fieldName,
	      fieldName,
	      vectorName);
      
      if(vectorType()==COMPLEX) {
	fprintf(outfile,"\t if(fscanf(infile,\"%%lf %%lf\",&_%s_%s[_i0].re,&_%s_%s[_i0].im) != 2) {\n",
		fieldName,vectorName,fieldName,vectorName);
      }
      else if(vectorType()==DOUBLE) {
	fprintf(outfile,"\t if(fscanf(infile,\"%%lf\",&_%s_%s[_i0]) != 1) {\n",fieldName,vectorName);
      }
      fprintf(outfile,
	      "\t\t printf(\"Error loading '%s' from file '%s': either bad float format or insufficient data\\n\");\n"
	      "\t\t _i0=_%s_size*_%s_%s_ncomponents;\n"
/*  This next line kills the program at this point.  
    We may or may not want to do this, although I just had a user 
    request that the program not continue past this error. (Wasted time)
    Reverting is as simple as deleting the next line.  (JJH)  */   
	      "\t\t exit(255);\n"
	      "\t }\n"
	      "\t _i0++;\n"
	      "}\n"
	      "\n"
	      "fclose(infile);\n",
	      vectorName,
	      myFileName.c_str(),
	      fieldName,
	      fieldName,
	      vectorName);
    }
    else if (myInputFileFormat == BINARY) {
      fprintf(outfile,
	      "FILE* infile = fopen(\"%s\",\"rb\");\n"
	      "\n"
	      "if (infile==0) {\n"
	      "\t printf(\"Error opening file '%s' for initialisation of vector '%s'\\n\");\n"
	      "\t return;\n"
	      "\t }\n"
	      "\n",
	      myFileName.c_str(),
	      myFileName.c_str(),
	      vectorName);
      
      if(vectorType()==COMPLEX) {
	fprintf(outfile,"\t complex _tempInput[_%s_size*_%s_%s_ncomponents];\n",
		fieldName,fieldName,vectorName);
	fprintf(outfile,"\t fread(&_tempInput, sizeof(complex), _%s_size*_%s_%s_ncomponents, infile);\n",
		fieldName,fieldName,vectorName);
	fprintf(outfile,"\t _%s_%s = _tempInput;\n",fieldName,vectorName);
      }
      else if(vectorType()==DOUBLE) {
	fprintf(outfile,"\t double _tempInput[_%s_size*_%s_%s_ncomponents];\n",
		fieldName,fieldName,vectorName);
	fprintf(outfile,"\t fread(&_tempInput, sizeof(double), _%s_size*_%s_%s_ncomponents, infile);\n",
		fieldName,fieldName,vectorName);
	fprintf(outfile,"\t _%s_%s = _tempInput;\n",fieldName,vectorName);
      }
      fprintf(outfile,
	      "\n"
	      "fclose(infile);\n");
    }
    else {
      // something must have really screwed up to get here...
      throw xmdsException("For some reason the input file is neither ascii nor binary...\nHow did we get here??");
    }
  }
  else {
    // initialisation from code
    if(simulation()->parameters()->stochastic) {
      fprintf(outfile,"const double _var = 1.0");
      for(unsigned long i=0;i<nDims;i++) {
	if(space(i)) {
	  fprintf(outfile,"/_%s_dk%li",fieldName,i);
	}
	else {
	  fprintf(outfile,"/_%s_dx%li",fieldName,i);
	}
      }
      fprintf(outfile,";\n");
      fprintf(outfile,"double *_noises = new double[_n_noises];\n");
      if(simulation()->parameters()->errorCheck) {
	fprintf(outfile,"double *_noises2 = new double[_n_noises];\n");
      }
      fprintf(outfile,"\n");
    }

    field()->vectors2space(outfile,initialSpace(),myVectorNamesList,"");

    list<XMLString> vectorNameList = myVectorNamesList;
    vectorNameList.push_back(*name());

    field()->openLoops(outfile,initialSpace(),vectorNameList);

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

    if(simulation()->parameters()->stochastic) {
      if(simulation()->parameters()->errorCheck) {
	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,"%sfor(unsigned long _s0=0;_s0<_n_noises;_s0++)\n",indent);
	fprintf(outfile,"%s	_noises[_s0] += _noises2[_s0];\n",indent);
	fprintf(outfile,"\n");
      }
      else {
	fprintf(outfile,"%s_make_noises(_gen,_var,_noises,_n_noises);\n",indent);
	fprintf(outfile,"\n");
      }
    }

    fprintf(outfile,"// ********** Code from vector element***********\n");
    fprintf(outfile,"%s\n",myCode.c_str());
    fprintf(outfile,"// **********************************************\n");
    fprintf(outfile,"\n");

    field()->closeLoops(outfile,initialSpace(),vectorNameList);
  }

  if(needsFFTWRoutines()) {
    fprintf(outfile,"\n");
    fprintf(outfile,"_%s_%s_space=%li;\n",field()->name()->c_str(),name()->c_str(),initialSpace());
  }
	
  if (simulation()->parameters()->stochastic) {
    fprintf(outfile, "    delete[] _noises;\n");
 	
    if (simulation()->parameters()->errorCheck) {
      fprintf(outfile, "    delete[] _noises2;\n");
    }
  }

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