xmdsintegrate.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: xmdsintegrate.cc,v 1.21 2004/10/06 07:52:38 joehope Exp $
*/

/*! @file xmdsintegrate.cc
  @brief Integrate element parsing classes and methods

  More detailed explanation...
*/

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

// ******************************************************************************
// ******************************************************************************
//                              xmdsIntegrate public
// ******************************************************************************
// ******************************************************************************

extern bool debugFlag;

long nxmdsIntegrates=0;  //!< Number of xmds integrate objects

// ******************************************************************************
xmdsIntegrate::xmdsIntegrate(
			     const xmdsSimulation *const yourSimulation,
			     const bool& yourVerboseMode) :
  xmdsSegment(yourSimulation,yourVerboseMode) {
  if(debugFlag) {
    nxmdsIntegrates++;
    printf("xmdsIntegrate::xmdsIntegrate\n");
    printf("nxmdsIntegrates=%li\n",nxmdsIntegrates);
  }
  myNextcoKey=0;
};

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

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


  list<XMLString> myXMLStringList;
  list<unsigned long> myULongList;

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

  // ************************************
  // find 'interval'

  getAssignmentStrings(yourElement,"interval",1,1,myXMLStringList);

  myInterval=*myXMLStringList.begin();

  if(verbose()) {
    printf("integrate interval = %s\n",myInterval.c_str());
  }

  // ************************************
  // find 'lattice'

  getAssignmentULongs(yourElement,"lattice",1,1,myULongList);

  myLattice=*myULongList.begin();

  if(myLattice==0) {
    throw xmdsException(yourElement,"Integration lattice must be >= 1 !");
  }

  if(verbose()) {
    printf("integrate lattice has %li points\n",myLattice);
  }

  // ************************************
  // find 'samples'

  getAssignmentULongs(yourElement,"samples",1,simulation()->output()->nMomentGroups(),mySamplesList);

  for(long unsigned int i=0;i<mySamplesList.size();i++) {
    if(samples(i) != 0) {
      if(myLattice != samples(i) * (myLattice/samples(i))) {
	throw xmdsException(yourElement,"moment group samples must be either 0 or else factors of the integration lattice");
      }
      if(verbose()) {
	printf("moment group #%li is sampled every %li points\n",i+1,myLattice/samples(i));
      }
    }
    else {
      if(verbose()) {
	printf("moment group #%li is not sampled\n",i+1);
      }
    }
  }

  // ************************************
  // find nonoises flag

  list<bool> tempNonoise;
  myNonoises=false;
  
  getAssignmentBools(yourElement,"no_noise",0,1,tempNonoise);
  
  if(tempNonoise.size()==1)
    {
    if(*tempNonoise.begin()) {
        myNonoises = true;
        if(verbose())
            printf("No noises will be definied in the filter.\n");
        }
  }

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

    // ************************************
    // find k_operators element

    const NodeList *const candidateElements = yourElement->getElementsByTagName("k_operators",0);

    if(candidateElements->length()>1) {
      throw xmdsException(yourElement,"Multiple <k_operators> elements defined");
    }

    if(candidateElements->length()==1) {

      // we have a k_operators element

      const Element *const kPropagationElement = dynamic_cast<const Element*> (candidateElements->item(0));

      if(verbose()) {
	printf("Processing <k_operators> element ...\n");
      }

      // ************************************
      // find constantK

      list<bool> aBoolList;

      getAssignmentBools(kPropagationElement,"constant",0,1,aBoolList);

      if(aBoolList.size()>0) {
	myConstantK=*aBoolList.begin();
	if(verbose()) {
	  if(myConstantK) {
	    printf("k-operators are propagation independent\n");
	  }
	  else {
	    printf("k-propagation are propagation dependent\n");
	  }
	}
      }
      else {
	myConstantK=1;
	printf("k-operators defaulting to propagation independent\n");
      }

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

      getAssignmentStrings(kPropagationElement,"vectors",0,0,myKVectorNamesList);

      if(myKVectorNamesList.size()==0) {
	// no vectors specified, therefore assume using only main vector
	myKVectorNamesList.push_back("main");
      }

      simulation()->field()->processVectors(myKVectorNamesList,simulation()->field()->geometry()->fullSpace());

      // ************************************
      // find operators

      getAssignmentStrings(kPropagationElement,"operator_names",1,0,myKOperatorNamesList);

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

      myKOperatorsCode=*kPropagationElement->textContent(0);

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

      if(myKOperatorsCode.isAllWhiteSpace()) {
	throw xmdsException(kPropagationElement,"No <k_operators> code defined!");
      }

      if(verbose()) {
	printf("<k_operators> code loaded\n");
      }
    }
  }

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

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

  // Vectors are now required by law in order to place the code element
  /*  if(myVectorNamesList.size()==0) {
  // no vectors specified, therefore assume using only main vector
  myVectorNamesList.push_back("main");
  }*/

  simulation()->field()->processVectors(myVectorNamesList,0);

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

  myPropagationCode=*yourElement->textContent(0);

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

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

  if(verbose()) {
    printf("processing propagation code ...\n");
  }


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

    // ************************************
    // find cross_propagation

    const NodeList *const candidateElements = yourElement->getElementsByTagName("cross_propagation",0);

    if(candidateElements->length() > 1) {
      throw xmdsException(yourElement,"only one <cross_propagation> allowed");
    }

    if(candidateElements->length() == 1) {

      // we have a cross_propagation element
                
      const Element *const crossPropagationElement = dynamic_cast<const Element*> (candidateElements->item(0));

      // ...but not if we're distributing the vector over many processors!
      if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
	throw xmdsException(crossPropagationElement,"I really can't handle using MPI when there are cross-propagation elements.");
      }

      // ************************************
      // find cross propagating vectors

      getAssignmentStrings(crossPropagationElement,"vectors",1,0,myCrossVectorNamesList);

      for(list<XMLString>::const_iterator pXMLString = myCrossVectorNamesList.begin(); pXMLString != myCrossVectorNamesList.end(); pXMLString++) {
	if(*pXMLString=="main") {
	  throw xmdsException(crossPropagationElement,"Cannot cross-propagate 'main' vector!");
	}
      }

      simulation()->field()->processVectors(myCrossVectorNamesList,0);

      // ************************************
      // find 'prop_dim'

      list<XMLString> anXMLStringList;

      getAssignmentStrings(crossPropagationElement,"prop_dim",1,1,anXMLStringList);

      // now find which member of the main field dimensions it is
      if(!simulation()->field()->geometry()->getDimNumber(*anXMLStringList.begin(),myCrossDimNumber)) {
	sprintf(errorMessage(),"prop_dim '%s' not a dimension of main field",anXMLStringList.begin()->c_str());
	throw xmdsException(crossPropagationElement,errorMessage());
      }

      if(verbose()) {
	printf("cross prop_dim is field dimension number %li\n",myCrossDimNumber);
      }

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

      myCrossPropagationCode=*crossPropagationElement->textContent(0);

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

      if(myCrossPropagationCode.isAllWhiteSpace()) {
	throw xmdsException(crossPropagationElement,"No <cross_propagation> code defined!");
      }

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

  // ************************************
  // process elements in order to process and place the code elements

  myNumNonLoopPropagation=0;
  myNumIntegrateMomentGroups=0;
  
  long unsigned int i;
  const NodeList* myElements;
  const Element* nextElement;

  myElements = yourElement->getElementsByTagName("*",0);
  i=0;

  while(i<myElements->length()) {
    nextElement = dynamic_cast<const Element*> (myElements->item(i));
    if(*nextElement->nodeName() == "functions") {
      XMLString someLoopPropagationCode = *nextElement->textContent(0);
      // check for non-white space charaters in code:
      if(someLoopPropagationCode.isAllWhiteSpace()) {
        throw xmdsException(nextElement,"No <functions> code defined!");
      }
      if(verbose()) {
        printf("    <functions> code loaded\n");
      }
      myNumNonLoopPropagation++;
      myNonLoopPropagationCodeList.push_back(someLoopPropagationCode);
      myCodeElementList.push_back(*nextElement->nodeName());
    }
    else if(*nextElement->nodeName() == "moment_group") {
      integrateMomentGroup tempIntegrateMG;
      
      XMLString someCode = *nextElement->textContent(0);
      // check for non-white space charaters in code:
      if(someCode.isAllWhiteSpace()) {
        throw xmdsException(nextElement,"No <moment_group> code defined!");
      }
      if(verbose()) {
        printf("   <moment_group> code loaded\n");
      }
      tempIntegrateMG.integrateMomentGroupCode = someCode;

      getAssignmentStrings(nextElement,"moments",1,0,tempIntegrateMG.momentNameList);      
      
      getAssignmentBools(nextElement,"integrate_dimension",1,simulation()->field()->geometry()->nDims(),tempIntegrateMG.integrateDimensionList);
      
      myNumIntegrateMomentGroups++;
      myIntegrateMomentGroupList.push_back(tempIntegrateMG);
      myCodeElementList.push_back(*nextElement->nodeName());
    }
    else if(*nextElement->nodeName() == "vectors") {
      myCodeElementList.push_back(*nextElement->nodeName());
    }

    i++;
  }
};

// ******************************************************************************
// ******************************************************************************
//                              xmdsIntegrate protected
// ******************************************************************************
// ******************************************************************************

// ******************************************************************************
void xmdsIntegrate::writeDefines(
				 FILE *const outfile) const {
  if(debugFlag) {
    printf("xmdsIntegrate::writeDefines\n");
  }

  if(usesKOperators()) {
    if(constantK()) {
      fprintf(outfile,"// integrate defines\n");
      fprintf(outfile,"\n");
      fprintf(outfile,"#define _segment%li_nkoperators %li\n",segmentNumber,(long)myKOperatorNamesList.size());
      fprintf(outfile,"\n");
    }
  }
};

// ******************************************************************************
void xmdsIntegrate::writeGlobals(
				 FILE *const outfile) const {
  if(debugFlag) {
    printf("xmdsIntegrate::writeGlobals\n");
  }

  if(usesKOperators()) {
    if(constantK()) {
      fprintf(outfile,"// integrate globals\n");
      fprintf(outfile,"\n");
      if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
	fprintf(outfile,"complex *_segment%li_k_operator_field;\n",segmentNumber);
      }
      else {
	fprintf(outfile,"complex *_segment%li_k_operator_field = new complex[_%s_size*_segment%li_nkoperators];\n",
		segmentNumber,simulation()->field()->name()->c_str(),segmentNumber);
      }
      fprintf(outfile,"\n");
    }
  }
};

// ******************************************************************************
void xmdsIntegrate::writePrototypes(
				    FILE *const outfile) const {
  if(debugFlag) {
    printf("xmdsIntegrate::writePrototypes\n");
  }

  writexSpacePrototype(outfile);
};

// ******************************************************************************
void xmdsIntegrate::writeRoutines(
				  FILE *const outfile) const {
  if(debugFlag) {
    printf("xmdsIntegrate::writeRoutines\n");
  }

  writexSpaceRoutine(outfile);
};

// ******************************************************************************
const XMLString* xmdsIntegrate::interval() const {
  if(debugFlag) {
    printf("xmdsIntegrate::interval\n");
  }

  return &myInterval;