/*
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.23 2005/05/02 02:35:58 sebwuester 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, const bool& adaptiveIP) :
xmdsSegment(yourSimulation,yourVerboseMode) {
myAdaptiveIP = adaptiveIP;
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 'smallmemory'
if(myAdaptiveIP){
list<bool> myBoolList;
getAssignmentBools(yourElement,"smallmemory",0,1,myBoolList);
if(myBoolList.size()==1) {
mySmallmemory=*myBoolList.begin();
if(verbose())
if(mySmallmemory) {
printf("\"save memory\" in k_propagate is 'yes'\n");
}
else {
printf("\"save memory\" in k _propagate 'no'\n");
}
}
else {
printf("\"save memory\" in k_propagate defaulting to no \n");
mySmallmemory=false;
}
}
// ************************************
// 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");
}
if(!myConstantK && myAdaptiveIP){
throw xmdsException("ARK45IP does not work with time-dependent k-operators. Use ARK45EX instead.");
}
// ************************************
// 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
// ******************************************************************************
// ******************************************************************************
// ******************************************************************************
const bool xmdsIntegrate::Smallmemory() const {
if(debugFlag) {
printf("xmdsIntegrate::Smallmemory\n");
}
return mySmallmemory;
};
// ******************************************************************************
void xmdsIntegrate::writeDefines(
FILE *const outfile) const {
if(debugFlag) {
printf("xmdsIntegrate::writeDefines\n");
}
if(usesKOperators()) {
if((myAdaptiveIP && !Smallmemory()) || (!myAdaptiveIP && 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(myAdaptiveIP){
if(usesKOperators())
if(!Smallmemory()) {
fprintf(outfile,"// integrate globals\n");
fprintf(outfile,"\n");
for(int n=1;n<6;n++)
if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
fprintf(outfile,"complex *_segment%li_k_operator_field_%i;\n",segmentNumber,n);
}
else {
fprintf(outfile,"complex *_segment%li_k_operator_field_%i = new complex[_%s_size*_segment%li_nkoperators];\n",
segmentNumber,n,simulation()->field()->name()->c_str(),segmentNumber);
}
fprintf(outfile,"\n");
}
}else{
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;
};
// ******************************************************************************
unsigned long xmdsIntegrate::lattice() const {
if(debugFlag) {
printf("xmdsIntegrate::lattice\n");
}
return myLattice;
};
// ******************************************************************************
unsigned long xmdsIntegrate::samples(
const unsigned long& index) const {
if(debugFlag) {
printf("xmdsIntegrate::KVectorNamesList\n");
}
if(index>=mySamplesList.size()) {
throw xmdsException("Internal range error in xmdsIntegrate::samples");
}
list<unsigned long>::const_iterator pULong = mySamplesList.begin();
for(unsigned long i=0; i<index; i++) {
pULong++;
}
return *pULong;
};
// ******************************************************************************
const list<XMLString>* xmdsIntegrate::KVectorNamesList() const {
if(debugFlag) {
printf("xmdsIntegrate::KVectorNamesList\n");
}
return &myKVectorNamesList;
};
// ******************************************************************************
bool xmdsIntegrate::constantK() const {
if(debugFlag) {
printf("xmdsIntegrate::constantK\n");
}
return myConstantK;
};
// ******************************************************************************
bool xmdsIntegrate::noNoises() const {
if(debugFlag) {
printf("xmdsIntegrate::noNoises\n");
}
return myNonoises;
};
// ******************************************************************************
unsigned long xmdsIntegrate::nKOperators() const {
if(debugFlag) {
printf("xmdsIntegrate::nKOperators\n");
}
return myKOperatorNamesList.size();
};
// ******************************************************************************
const XMLString* xmdsIntegrate::KOperator(
const unsigned long& index) const {
if(debugFlag) {
printf("xmdsIntegrate::KOperator\n");
}
if(index>=myKOperatorNamesList.size()) {
throw xmdsException("Internal range error in xmdsIntegrate::KOperator");
}
list<XMLString>::const_iterator pXMLString = myKOperatorNamesList.begin();
for(unsigned long i=0; i<index; i++) {
pXMLString++;
}
return &*pXMLString;
};
// ******************************************************************************
bool xmdsIntegrate::getKOperator(
const XMLString& operatorName,
unsigned long& index) const {
if(debugFlag) {
printf("xmdsField::getVector\n");
}
index=0;
for(list<XMLString>::const_iterator pXMLString = myKOperatorNamesList.begin(); pXMLString != myKOperatorNamesList.end(); pXMLString++) {
if(*pXMLString==operatorName) {
return 1;
}
index++;
}
return 0;
};
// ******************************************************************************
const XMLString* xmdsIntegrate::KOperatorsCode() const {
if(debugFlag) {
printf("xmdsIntegrate::KOperatorsCode\n");
}
return &myKOperatorsCode;
};
// ******************************************************************************
list<XMLString>* xmdsIntegrate::vectorNamesList() {
if(debugFlag) {
printf("xmdsIntegrate::vectorNamesList\n");
}
return &myVectorNamesList;
};
// ******************************************************************************
const list<XMLString>* xmdsIntegrate::vectorNamesList() const {
if(debugFlag) {
printf("xmdsIntegrate::vectorNamesList\n");
}
return &myVectorNamesList;
};
// ******************************************************************************
XMLString* xmdsIntegrate::propagationCode() {
if(debugFlag) {
printf("xmdsIntegrate::propagationCode\n");
}
return &myPropagationCode;
};
// ******************************************************************************
const XMLString* xmdsIntegrate::propagationCode() const {
if(debugFlag) {
printf("xmdsIntegrate::propagationCode\n");
}
return &myPropagationCode;
};
// ******************************************************************************
const list<XMLString>* xmdsIntegrate::crossVectorNamesList() const {
if(debugFlag) {
printf("xmdsIntegrate::crossVectorNamesList\n");
}
return &myCrossVectorNamesList;
};
// ******************************************************************************
unsigned long xmdsIntegrate::crossDimNumber() const {
if(debugFlag) {
printf("xmdsIntegrate::crossDimNumber\n");
}
return myCrossDimNumber;
};
// ******************************************************************************
const XMLString* xmdsIntegrate::crossPropagationCode() const {
if(debugFlag) {
printf("xmdsIntegrate::crossPropagationCode\n");
}
return &myCrossPropagationCode;
};
// ******************************************************************************
list<XMLString>* xmdsIntegrate::nonLoopPropagationCodeList() {
if(debugFlag) {
printf("xmdsIntegrate::nonLoopPropagationCodeList\n");
}
return &myNonLoopPropagationCodeList;
};
// ******************************************************************************
const list<XMLString>* xmdsIntegrate::nonLoopPropagationCodeList() const {
if(debugFlag) {
printf("xmdsIntegrate::nonLoopPropagationCodeList\n");
}
return &myNonLoopPropagationCodeList;
};
// ******************************************************************************
list<XMLString>* xmdsIntegrate::codeElementList() {
if(debugFlag) {
printf("xmdsIntegrate::codeElementList\n");
}
return &myCodeElementList;
};
// ******************************************************************************
const list<XMLString>* xmdsIntegrate::codeElementList() const {
if(debugFlag) {
printf("xmdsIntegrate::codeElementList\n");
}
return &myCodeElementList;
};
// ******************************************************************************
bool xmdsIntegrate::usesKOperators() const {
if(debugFlag) {
printf("xmdsIntegrate::usesKOperators\n");
}
return (mycoStructList.size()>0);
};
// ******************************************************************************
const bool xmdsIntegrate::AdaptiveIP() const {
if(debugFlag) {
printf("xmdsIntegrate::AdaptiveIP\n");
}
return myAdaptiveIP;
};
// ******************************************************************************
long unsigned int xmdsIntegrate::numNonLoopPropagation() const {
if(debugFlag) {
printf("xmdsIntegrate::numNonLoopPropagation\n");
}
return (myNumNonLoopPropagation);
};
// ******************************************************************************
long unsigned int xmdsIntegrate::numIntegrateMomentGroups() const {
if(debugFlag) {
printf("xmdsIntegrate::numIntegrateMomentGroups\n");
}
return (myNumIntegrateMomentGroups);
};
// ******************************************************************************
list<integrateMomentGroup>* xmdsIntegrate::integrateMomentGroupList() {
if(debugFlag) {
printf("xmdsIntegrate::integrateMomentGroupList\n");
}
return &myIntegrateMomentGroupList;
};
// ******************************************************************************
const list<integrateMomentGroup>* xmdsIntegrate::integrateMomentGroupList() const {
if(debugFlag) {
printf("xmdsIntegrate::integrateMomentGroupList\n");
}
return &myIntegrateMomentGroupList;
};
// ******************************************************************************
bool xmdsIntegrate::getcoStruct(
const unsigned long& componentNumber,
const coStruct*& thecoStruct) const {
if(debugFlag) {
printf("xmdsIntegrate::getcoStruct\n");
}
for(list<coStruct>::const_iterator pcoStruct = mycoStructList.begin(); pcoStruct != mycoStructList.end(); pcoStruct++) {
if(pcoStruct->componentNumber == componentNumber) {
thecoStruct = &*pcoStruct;
return 1;
}
}
return 0;
};
// ******************************************************************************
bool xmdsIntegrate::getcoKey(
const unsigned long& componentNumber,
const unsigned long& operatorNumber,
unsigned long& coKey) const {
if(debugFlag) {
printf("xmdsIntegrate::getcoPairKey\n");
}
const coStruct* thecoStruct;
if(!getcoStruct(componentNumber,thecoStruct)) {
return 0;
}
list<unsigned long>::const_iterator pULong1 = thecoStruct->operatorNumbersList.begin();
list<unsigned long>::const_iterator pULong2 = thecoStruct->coKeysList.begin();
while(pULong1 != thecoStruct->operatorNumbersList.end()) {
if(*pULong1==operatorNumber) {
coKey = *pULong2;
return 1;
}
pULong1++;
pULong2++;
}
return 0;
};
// ******************************************************************************
unsigned long xmdsIntegrate::addcoPair(
const unsigned long& componentNumber,
const unsigned long& operatorNumber) {
if(debugFlag) {
printf("xmdsIntegrate::addcoPair\n");
}
coStruct* existingcoStruct = 0;
for(list<coStruct>::iterator pcoStruct = mycoStructList.begin(); pcoStruct != mycoStructList.end(); pcoStruct++) {
if(pcoStruct->componentNumber == componentNumber) {
existingcoStruct = &*pcoStruct;
}
}
if(existingcoStruct != 0) {
existingcoStruct->operatorNumbersList.push_back(operatorNumber);
existingcoStruct->coKeysList.push_back(myNextcoKey);
}
else {
coStruct nextcoStruct;
nextcoStruct.componentNumber = componentNumber;
nextcoStruct.operatorNumbersList.push_back(operatorNumber);
nextcoStruct.coKeysList.push_back(myNextcoKey);
mycoStructList.push_back(nextcoStruct);
}
myNextcoKey++;
return myNextcoKey-1;
};
// ******************************************************************************
bool xmdsIntegrate::findNextcoPair(
XMLString& operatorName,
XMLString& componentName,
unsigned long& start,
unsigned long& end) const {
if(debugFlag) {
printf("xmdsIntegrate::findNextcoPair\n");
}
// here we are looking for the next **** operatorName S? '[' S? componentName S? ']' **** production
// starting from the offset provided
const char* code=myPropagationCode.c_str();
end = start;
unsigned long openBrace;
unsigned long nameStart;
unsigned long nameEnd;
// find next [
while((code[end] != '[') && (code[end] != 0)) {
end++;
}
if(code[end] == 0) {
return 0;
}
// ok, we found the first [
openBrace = end;
// now skip white space
end++;
while(XMLChar::isWhiteSpace(code[end]) && (code[end] != 0)) {
end++;
}
if(code[end] == 0) {
return 0;
}
if(!(XMLChar::isLatinLetter(code[end]) | (code[end]=='_'))) {
return 0;
}
nameStart = end;
end++;
while((XMLChar::isLatinLetter(code[end]) | XMLChar::isLatinDigit(code[end]) | (code[end]=='_'))&(code[end] != 0)) {
end++;
}
if(code[end] == 0) {
return 0;
}
nameEnd = end;
//check that we have a ] on the end
// skip white space
while(XMLChar::isWhiteSpace(code[end]) && (code[end] != 0)) {
end++;
}
if(code[end] != ']') {
return 0;
}
// ok, we definititely had a nice clean [name]
end++;
myPropagationCode.subString(componentName,nameStart,nameEnd);
// now work backwards to find the operator name
start=openBrace-1;
// skip any white space
while(XMLChar::isWhiteSpace(code[start])&(start>0)) {
start--;
}
if(!(XMLChar::isLatinLetter(code[start]) | XMLChar::isLatinDigit(code[start]) | (code[start]=='_'))) {
return 0;
}
nameEnd=start+1;
while((XMLChar::isLatinLetter(code[start]) | XMLChar::isLatinDigit(code[start]) | (code[start]=='_'))&(start>0)) {
start--;
}
if(!(XMLChar::isLatinLetter(code[start]) | XMLChar::isLatinDigit(code[start]) | (code[start]=='_'))) {
start++;
}
if(!(XMLChar::isLatinLetter(code[start]) | (code[start]=='_'))) {
return 0;
}
nameStart=start;
myPropagationCode.subString(operatorName,nameStart,nameEnd);
return 1;
};
// ******************************************************************************
void xmdsIntegrate::writexSpacePrototype(
FILE *const outfile) const {
if(debugFlag) {
printf("xmdsIntegrate::writexSpacePrototype\n");
}
fprintf(outfile,"void _segment%li_calculate_delta_a(\n",segmentNumber);
fprintf(outfile," const double& _step");
fprintf(outfile,",\n const unsigned long cycle");
if((simulation()->parameters()->stochastic)&&(!myNonoises)) {
fprintf(outfile,",\n const double *const _noise_vector");
}
fprintf(outfile,");\n");
fprintf(outfile,"\n");
};
// ******************************************************************************
void xmdsIntegrate::writexSpaceRoutine(
FILE *const outfile) const {
if(debugFlag) {
printf("xmdsIntegrate::writexSpaceRoutine\n");
}
const unsigned long nDims = simulation()->field()->geometry()->nDims();
const char *const fieldName = simulation()->field()->name()->c_str();
fprintf(outfile,"// *************************\n");
fprintf(outfile,"void _segment%li_calculate_delta_a(\n",segmentNumber);
fprintf(outfile," const double& _step");
fprintf(outfile,",\n const unsigned long cycle");
if((simulation()->parameters()->stochastic)&&(!myNonoises)) {
fprintf(outfile,",\n const double *const _noise_vector");
}
fprintf(outfile,") {\n");
fprintf(outfile,"\n");
const xmdsVector* mainVector;
if(!simulation()->field()->getVector("main",mainVector)) {
throw xmdsException("Internal error in xmdsIntegrate::writexSpaceRoutine: cannot find 'main' vector");
}
const char* typeName="";
if(mainVector->vectorType()==COMPLEX) {
typeName="complex";
}
else if(mainVector->vectorType()==DOUBLE) {
typeName="double";
}
for(unsigned long i=0;i<mainVector->nComponents();i++) {
fprintf(outfile,"%s d%s_d%s;\n",typeName,mainVector->componentName(i)->c_str(),simulation()->parameters()->propDimName.c_str());
}
fprintf(outfile,"\n");
// add cross vectors to total vectors to use, but no duplicating names!
list<XMLString> myTotalVectorsList = *vectorNamesList();
for(list<XMLString>::const_iterator pXMLString = crossVectorNamesList()->begin(); pXMLString != crossVectorNamesList()->end(); pXMLString++) {
list<XMLString>::const_iterator pXMLString2 = myTotalVectorsList.begin();
while((pXMLString2 != myTotalVectorsList.end()) && (*pXMLString2 != *pXMLString)) {
pXMLString2++;
}
if(*pXMLString2 != *pXMLString) {
myTotalVectorsList.push_back(*pXMLString);
}
}
simulation()->field()->vectors2space(outfile,0,myTotalVectorsList,"");
if(crossVectorNamesList()->size() > 0) {
fprintf(outfile,"_segment%li_calculate_cross_field();\n",segmentNumber);
fprintf(outfile,"\n");
}
if((simulation()->parameters()->stochastic)&&(!myNonoises)) {
fprintf(outfile,"const double* _noises = _noise_vector;\n");
fprintf(outfile,"\n");
}
list<XMLString>::const_iterator nextNLPElement = myNonLoopPropagationCodeList.begin();
list<integrateMomentGroup>::const_iterator nextMGElement = myIntegrateMomentGroupList.begin();
list<XMLString> theCodeList = *codeElementList();
if(!(theCodeList.size()==1+numIntegrateMomentGroups()+numNonLoopPropagation())) {
throw xmdsException("The list of integrate code elements is the wrong length!");
}
long unsigned int whichMG = 0;
list<long> deleteMGArrayList;
for(list<XMLString>::const_iterator codeElement = theCodeList.begin(); codeElement != theCodeList.end(); codeElement++) {
if(!strcmp(codeElement->c_str(),"vectors")) {
simulation()->field()->openLoops(outfile,0,myTotalVectorsList);
fprintf(outfile,"\n");
char indent[64];
for(unsigned long i=0;i<nDims;i++) {
indent[i]=0x09;
}
indent[nDims]=0;
// integrate moment group pointer definition
long unsigned int tempWhichMG = 0;
for(list<integrateMomentGroup>::const_iterator tempNextMGElement = myIntegrateMomentGroupList.begin();
tempNextMGElement != myIntegrateMomentGroupList.end(); tempNextMGElement++) {
long unsigned int nDimsRemaining = simulation()->field()->geometry()->nDims();
tempWhichMG++;
for(list<bool>::const_iterator nextIntegrateDimension = tempNextMGElement->integrateDimensionList.begin();
nextIntegrateDimension != tempNextMGElement->integrateDimensionList.end(); nextIntegrateDimension++){
if(*nextIntegrateDimension) {
nDimsRemaining--;
}
}
if(nDimsRemaining>0) {
fprintf(outfile," long _img%li_pointer = ",tempWhichMG);
for(long unsigned int j=0; j<nDimsRemaining-1; j++) {
fprintf(outfile,"(");
}
long unsigned int whichMoment = 0;
list<bool>::const_iterator nextIntegrateDimension = tempNextMGElement->integrateDimensionList.begin();
while(*nextIntegrateDimension) {
nextIntegrateDimension++;
whichMoment++;
}
fprintf(outfile,"_i%li",whichMoment);
nextIntegrateDimension++;
whichMoment++;
for(long unsigned int k=0; k<nDimsRemaining-1; k++) {
while(*nextIntegrateDimension) {
nextIntegrateDimension++;
whichMoment++;
}
fprintf(outfile,"*_%s_lattice%li+_i%li)",fieldName,whichMoment,whichMoment);
nextIntegrateDimension++;
whichMoment++;
}
fprintf(outfile,"*%i;\n",tempNextMGElement->momentNameList.size());
}
}
fprintf(outfile,"// ************** propagation code **************\n");
fprintf(outfile,"%s\n",propagationCode()->c_str());
fprintf(outfile,"// **********************************************\n");
fprintf(outfile,"\n");
for(long unsigned int i=0;i<mainVector->nComponents();i++) {
fprintf(outfile,"%s_active_%s_main[_%s_main_index_pointer + %li] = d%s_d%s*_step;\n",
indent,fieldName,fieldName,i,mainVector->componentName(i)->c_str(),simulation()->parameters()->propDimName.c_str());
}
fprintf(outfile,"\n");
if((simulation()->parameters()->stochastic)&&(!myNonoises)) {
fprintf(outfile,"%s_noises += _n_noises;\n",indent);
fprintf(outfile,"\n");
}
simulation()->field()->closeLoops(outfile,0,myTotalVectorsList);
fprintf(outfile,"\n");
}
else if(!strcmp(codeElement->c_str(),"functions")) {
fprintf(outfile,"// ************** propagation code **************\n");
fprintf(outfile,"%s\n",nextNLPElement->c_str());
fprintf(outfile,"// **********************************************\n");
fprintf(outfile,"\n");
nextNLPElement++;
}
else if(!strcmp(codeElement->c_str(),"moment_group")) {
whichMG++;
long unsigned int nDimsRemaining = simulation()->field()->geometry()->nDims();
for(list<bool>::const_iterator nextIntegrateDimension = nextMGElement->integrateDimensionList.begin();
nextIntegrateDimension != nextMGElement->integrateDimensionList.end(); nextIntegrateDimension++){
if(*nextIntegrateDimension) {
nDimsRemaining--;
}
}
if(nDimsRemaining>0) {
deleteMGArrayList.push_back(whichMG);
}
fprintf(outfile,"// ************** integrate moment group code **************\n");
//setup defines for the integrate_moment_groups if they are arrays
long unsigned int whichMoment=0;
if(nDimsRemaining>0) {
for(list<XMLString>::const_iterator nextMomentName = nextMGElement->momentNameList.begin();
nextMomentName != nextMGElement->momentNameList.end(); nextMomentName++){
fprintf(outfile,"#define %s _img%li_array[_img%li_pointer+%li]\n",nextMomentName->c_str(),whichMG,whichMG,whichMoment);
whichMoment++;
}
}
fprintf(outfile,"\n");
//setup actual arrays for the integrate_moment_groups (MPI and non-MPI)
if(nDimsRemaining>0) {
fprintf(outfile,"complex *_img%li_array = new complex[%i",whichMG,nextMGElement->momentNameList.size());
if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
whichMoment=0;
for(list<bool>::const_iterator nextIntegrateDimension = nextMGElement->integrateDimensionList.begin();
nextIntegrateDimension != nextMGElement->integrateDimensionList.end(); nextIntegrateDimension++) {
if(!*nextIntegrateDimension) {
if(whichMoment==0){
fprintf(outfile,"*local_nx");
}
else {
fprintf(outfile,"*_%s_lattice%li",fieldName,whichMoment);
}
}
whichMoment++;
}
}
else {
whichMoment=0;
for(list<bool>::const_iterator nextIntegrateDimension = nextMGElement->integrateDimensionList.begin();
nextIntegrateDimension != nextMGElement->integrateDimensionList.end(); nextIntegrateDimension++) {
if(!*nextIntegrateDimension) {
fprintf(outfile,"*_%s_lattice%li",fieldName,whichMoment);
}
whichMoment++;
}
}
fprintf(outfile,"];\n");
}
else {
for(list<XMLString>::const_iterator nextMomentName = nextMGElement->momentNameList.begin();
nextMomentName != nextMGElement->momentNameList.end(); nextMomentName++) {
fprintf(outfile,"complex %s = rcomplex(0,0);\n",nextMomentName->c_str());
}
}
fprintf(outfile,"\n");
//zero the arrays and setup a local environment for the sum
if(nDimsRemaining>0) {
fprintf(outfile,"for(long _i0=0; _i0<%i",nextMGElement->momentNameList.size());
if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
whichMoment=0;
for(list<bool>::const_iterator nextIntegrateDimension = nextMGElement->integrateDimensionList.begin();
nextIntegrateDimension != nextMGElement->integrateDimensionList.end(); nextIntegrateDimension++) {
if(!*nextIntegrateDimension) {
if(whichMoment==0){
fprintf(outfile,"*local_nx");
}
else {
fprintf(outfile,"*_%s_lattice%li",fieldName,whichMoment);
}
}
whichMoment++;
}
}
else {
whichMoment=0;
for(list<bool>::const_iterator nextIntegrateDimension = nextMGElement->integrateDimensionList.begin();
nextIntegrateDimension != nextMGElement->integrateDimensionList.end(); nextIntegrateDimension++) {
if(!*nextIntegrateDimension) {
fprintf(outfile,"*_%s_lattice%li",fieldName,whichMoment);
}
whichMoment++;
}
}
fprintf(outfile,"; _i0++) {\n");
fprintf(outfile," _img%li_array[_i0] = rcomplex(0,0);\n",whichMG);
fprintf(outfile,"}\n");
}
fprintf(outfile,"{\n");
// open loops
simulation()->field()->openLoops(outfile,0,myTotalVectorsList);
// pointer definition
if(nDimsRemaining>0) {
fprintf(outfile," long _img%li_pointer = ",whichMG);
for(long unsigned int j=0; j<nDimsRemaining-1; j++) {
fprintf(outfile,"(");
}
whichMoment = 0;
list<bool>::const_iterator nextIntegrateDimension = nextMGElement->integrateDimensionList.begin();
while(*nextIntegrateDimension) {
nextIntegrateDimension++;
whichMoment++;
}
fprintf(outfile,"_i%li",whichMoment);
nextIntegrateDimension++;
whichMoment++;
for(long unsigned int k=0; k<nDimsRemaining-1; k++) {
while(*nextIntegrateDimension) {
nextIntegrateDimension++;
whichMoment++;
}
fprintf(outfile,"*_%s_lattice%li+_i%li)",fieldName,whichMoment,whichMoment);
nextIntegrateDimension++;
whichMoment++;
}
fprintf(outfile,"*%i;\n",nextMGElement->momentNameList.size());
}
// code
fprintf(outfile,"%s\n",nextMGElement->integrateMomentGroupCode.c_str());
// close loops
simulation()->field()->closeLoops(outfile,0,myTotalVectorsList);
fprintf(outfile,"}\n");
// if there was any integration, multiply by the relevant volume element
if(nDimsRemaining<simulation()->field()->geometry()->nDims()) {
if(nDimsRemaining>0) {
fprintf(outfile,"for(long _i0=0; _i0<%i",nextMGElement->momentNameList.size());
if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
whichMoment=0;
for(list<bool>::const_iterator nextIntegrateDimension = nextMGElement->integrateDimensionList.begin();
nextIntegrateDimension != nextMGElement->integrateDimensionList.end(); nextIntegrateDimension++) {
if(!*nextIntegrateDimension) {
if(whichMoment==0){
fprintf(outfile,"*local_nx");
}
else {
fprintf(outfile,"*_%s_lattice%li",fieldName,whichMoment);
}
}
whichMoment++;
}
}
else {
whichMoment=0;
for(list<bool>::const_iterator nextIntegrateDimension = nextMGElement->integrateDimensionList.begin();
nextIntegrateDimension != nextMGElement->integrateDimensionList.end(); nextIntegrateDimension++) {
if(!*nextIntegrateDimension) {
fprintf(outfile,"*_%s_lattice%li",fieldName,whichMoment);
}
whichMoment++;
}
}
whichMoment=0;
fprintf(outfile,"; _i0++) {\n");
fprintf(outfile," _img%li_array[_i0] = _img%li_array[_i0]",whichMG,whichMG);
for(list<bool>::const_iterator nextIntegrateDimension = nextMGElement->integrateDimensionList.begin();
nextIntegrateDimension != nextMGElement->integrateDimensionList.end(); nextIntegrateDimension++){
if(*nextIntegrateDimension) {
fprintf(outfile,"*_%s_dx%li",fieldName,whichMoment);
}
whichMoment++;
}
fprintf(outfile,";\n");
fprintf(outfile," }\n");
}
else {
for(list<XMLString>::const_iterator nextMomentName = nextMGElement->momentNameList.begin();
nextMomentName != nextMGElement->momentNameList.end(); nextMomentName++) {
whichMoment = 0;
fprintf(outfile," %s = %s",nextMomentName->c_str(),nextMomentName->c_str());
for(list<bool>::const_iterator nextIntegrateDimension = nextMGElement->integrateDimensionList.begin();
nextIntegrateDimension != nextMGElement->integrateDimensionList.end(); nextIntegrateDimension++){
if(*nextIntegrateDimension) {
fprintf(outfile,"*_%s_dx%li",fieldName,whichMoment);
}
whichMoment++;
}
fprintf(outfile,";\n");
}
}
}
// if mpi and first variable integrated then mpireduce
if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
if(* nextMGElement->integrateDimensionList.begin()) {
if(nDimsRemaining>0) {
fprintf(outfile,"\n complex *_temp_img%li_array = new complex[%i",whichMG,nextMGElement->momentNameList.size());
whichMoment=0;
for(list<bool>::const_iterator nextIntegrateDimension = nextMGElement->integrateDimensionList.begin();
nextIntegrateDimension != nextMGElement->integrateDimensionList.end(); nextIntegrateDimension++) {
if(!*nextIntegrateDimension) {
fprintf(outfile,"*_%s_lattice%li",fieldName,whichMoment);
}
whichMoment++;
}
fprintf(outfile,"];\n");
fprintf(outfile," MPI_Reduce(_img%li_array,_temp_img%li_array,2*%i",whichMG,whichMG,nextMGElement->momentNameList.size());
whichMoment=0;
for(list<bool>::const_iterator nextIntegrateDimension = nextMGElement->integrateDimensionList.begin();
nextIntegrateDimension != nextMGElement->integrateDimensionList.end(); nextIntegrateDimension++) {
if(!*nextIntegrateDimension) {
fprintf(outfile,"*_%s_lattice%li",fieldName,whichMoment);
}
whichMoment++;
}
fprintf(outfile,",MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);\n");
fprintf(outfile,"\n delete[] _temp_img%li_array;\n",whichMG);
}
else {
fprintf(outfile,"\n complex *_temp%li_mpireduce = new complex;\n",whichMG);
for(list<XMLString>::const_iterator nextMomentName = nextMGElement->momentNameList.begin();
nextMomentName != nextMGElement->momentNameList.end(); nextMomentName++) {
fprintf(outfile," MPI_Reduce(&%s,_temp%li_mpireduce,2,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);\n",nextMomentName->c_str(),whichMG);
fprintf(outfile,"\n delete _temp%li_mpireduce;\n",whichMG);
}
}
}
}
fprintf(outfile,"\n");
fprintf(outfile,"// *********************************************************\n");
fprintf(outfile,"\n");
nextMGElement++;
}
else {
throw xmdsException("Unknown code element in the integrate block!");
}
}
// delete any arrays that were generated
for(list<long>::const_iterator deleteMGArray = deleteMGArrayList.begin();
deleteMGArray != deleteMGArrayList.end(); deleteMGArray++) {
fprintf(outfile,"delete[] _img%li_array;\n",*deleteMGArray);
}
fprintf(outfile,"}\n");
fprintf(outfile,"\n");
};
// ******************************************************************************
// ******************************************************************************
// xmdsIntegrate private
// ******************************************************************************
// ******************************************************************************
// ******************************************************************************
void xmdsIntegrate::outputSampleCount() const {
if(debugFlag) {
printf("xmdsIntegrate::outputSampleCount\n");
}
simulation()->output()->addSamples(mySamplesList);
};
// ******************************************************************************
void xmdsIntegrate::writeInitialisationCalls(
FILE *const outfile) const {
if(debugFlag) {
printf("xmdsIntegrate::writeInitialisationCalls\n");
}
if(usesKOperators()&&!Smallmemory()) {
if(myAdaptiveIP){
for(int i=1;i<6;i++)
if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
fprintf(outfile,"_segment%li_k_operator_field_%i = new complex[total_local_size*_segment%li_nkoperators];\n",
segmentNumber,i,segmentNumber);
fprintf(outfile,"\n");
}
}else{
if(constantK())
fprintf(outfile,"_segment%li_calculate_k_operator_field();\n",segmentNumber);
}
}
};
syntax highlighted by Code2HTML, v. 0.9.1