/*
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: xmdsintegratesiip.cc,v 1.25 2005/07/27 06:57:43 joehope Exp $
*/
/*! @file xmdsintegratesiip.cc
@brief Integrate element parsing classes and methods; semi-implicit method in the interaction picture
More detailed explanation...
*/
#include<xmlbasics.h>
#include<dom3.h>
#include<xmdsutils.h>
#include<xmdsclasses.h>
// ******************************************************************************
// ******************************************************************************
// xmdsIntegrateSIIP public
// ******************************************************************************
// ******************************************************************************
extern bool debugFlag;
long nxmdsIntegrateSIIPs=0; //!< Number of xmds integrate SIIP objects
// ******************************************************************************
xmdsIntegrateSIIP::xmdsIntegrateSIIP(
const xmdsSimulation *const yourSimulation,
const bool& yourVerboseMode) :
xmdsIntegrate(yourSimulation,yourVerboseMode),
xmdsIntegrateIP(yourSimulation,yourVerboseMode),
xmdsIntegrateSI(yourSimulation,yourVerboseMode) {
if(debugFlag) {
nxmdsIntegrateSIIPs++;
printf("xmdsIntegrateSIIP::xmdsIntegrateSIIP\n");
printf("nxmdsIntegrateSIIPs=%li\n",nxmdsIntegrateSIIPs);
}
};
// ******************************************************************************
xmdsIntegrateSIIP::~xmdsIntegrateSIIP() {
if(debugFlag) {
nxmdsIntegrateSIIPs--;
printf("xmdsIntegrateSIIP::~xmdsIntegrateSIIP\n");
printf("nxmdsIntegrateSIIPs=%li\n",nxmdsIntegrateSIIPs);
}
};
// ******************************************************************************
void xmdsIntegrateSIIP::processElement(
const Element *const yourElement) {
if(debugFlag) {
printf("xmdsIntegrateSIIP::processElement\n");
}
if(verbose()) {
printf("Processing integrate SIIP element ...\n");
}
xmdsIntegrate::processElement(yourElement);
xmdsIntegrateIP::processElement(yourElement);
xmdsIntegrateSI::processElement(yourElement);
};
// ******************************************************************************
// ******************************************************************************
// xmdsIntegrateSIIP private
// ******************************************************************************
// ******************************************************************************
// ******************************************************************************
void xmdsIntegrateSIIP::writePrototypes(
FILE *const outfile) const {
if(debugFlag) {
printf("xmdsIntegrateSIIP::writePrototypes\n");
}
fprintf(outfile,"// ********************************************************\n");
fprintf(outfile,"// segment %li (SIIP) prototypes\n",segmentNumber);
fprintf(outfile,"\n");
xmdsIntegrate::writePrototypes(outfile);
xmdsIntegrateIP::writePrototypes(outfile);
fprintf(outfile,"void _segment%li(unsigned long cycle);\n",segmentNumber);
fprintf(outfile,"\n");
};
// ******************************************************************************
void xmdsIntegrateSIIP::writexSpacePrototype(
FILE *const outfile) const {
if(debugFlag) {
printf("xmdsIntegrateSIIP::writexSpacePrototype\n");
}
fprintf(outfile,"void _segment%li_x_propagate(\n",segmentNumber);
fprintf(outfile," const double& _step");
fprintf(outfile,",\n const unsigned long cycle");
if((simulation()->parameters()->stochastic)&&(!noNoises())) {
fprintf(outfile,",\n const unsigned long& _generator_flag");
}
fprintf(outfile,");\n");
fprintf(outfile,"\n");
};
// ******************************************************************************
void xmdsIntegrateSIIP::writeRoutines(
FILE *const outfile) const {
if(debugFlag) {
printf("xmdsIntegrateSIIP::writeRoutines\n");
}
fprintf(outfile,"// ********************************************************\n");
fprintf(outfile,"// segment %li (SIIP) routines\n",segmentNumber);
fprintf(outfile,"\n");
xmdsIntegrate::writeRoutines(outfile);
xmdsIntegrateIP::writeRoutines(outfile);
writeMainIntegrateRoutine(outfile);
};
// ******************************************************************************
void xmdsIntegrateSIIP::writeMainIntegrateRoutine(
FILE *const outfile) const {
if(debugFlag) {
printf("xmdsIntegrateSIIP::writeMainIntegrateRoutine\n");
}
fprintf(outfile,"/* ******************************************** */\n");
fprintf(outfile,"void _segment%li(unsigned long cycle) {\n",segmentNumber);
fprintf(outfile,"\n");
fprintf(outfile,"for(unsigned long _i0=0;_i0<%li;_i0++) {\n",lattice());
fprintf(outfile,"\n");
if(simulation()->parameters()->errorCheck) {
fprintf(outfile," if(_half_step) {\n");
fprintf(outfile,"\n");
fprintf(outfile," const double _step = 0.5*%s/(double)%li;\n",interval()->c_str(),lattice());
fprintf(outfile,"\n");
writeSingleStepCode(outfile,FIRST_HALFSTEP);
fprintf(outfile,"\n");
writeSingleStepCode(outfile,SECOND_HALFSTEP);
fprintf(outfile," }\n");
fprintf(outfile," else {\n");
fprintf(outfile,"\n");
fprintf(outfile," const double _step = %s/(double)%li;\n",interval()->c_str(),lattice());
fprintf(outfile,"\n");
writeSingleStepCode(outfile,FULLSTEP);
fprintf(outfile," }\n");
}
else {
fprintf(outfile,"\n");
fprintf(outfile," const double _step = %s/(double)%li;\n",interval()->c_str(),lattice());
fprintf(outfile,"\n");
writeSingleStepCode(outfile,FULLSTEP);
}
for(long unsigned int i=0;i<simulation()->output()->nMomentGroups();i++) {
if(samples(i)!=0){
fprintf(outfile,"\n");
fprintf(outfile," if(%li*((_i0+1)/%li)==(_i0+1))\n",lattice()/samples(i),lattice()/samples(i));
fprintf(outfile," _mg%li_sample();\n",i);
}
}
fprintf(outfile," }\n");
fprintf(outfile,"}\n");
fprintf(outfile,"\n");
};
// ******************************************************************************
void xmdsIntegrateSIIP::writexSpaceRoutine(
FILE *const outfile) const {
if(debugFlag) {
printf("xmdsIntegrateSIIP::writexSpaceRoutine\n");
}
const long unsigned int nDims = simulation()->field()->geometry()->nDims();
const char *const fieldName = simulation()->field()->name()->c_str();
fprintf(outfile,"// *************************\n");
fprintf(outfile,"void _segment%li_x_propagate(\n",segmentNumber);
fprintf(outfile," const double& _step");
fprintf(outfile,",\n const unsigned long cycle");
if((simulation()->parameters()->stochastic)&&(!noNoises())) {
fprintf(outfile,",\n const unsigned long& _generator_flag");
}
fprintf(outfile,") {\n");
fprintf(outfile,"\n");
const xmdsVector* mainVector;
if(!simulation()->field()->getVector("main",mainVector)) {
throw xmdsException("Internal error in xmdsIntegrateSIIP::writexSpaceRoutine: cannot find 'main' vector");
}
const char* typeName="";
if(mainVector->vectorType()==COMPLEX) {
typeName="complex";
}
else if(mainVector->vectorType()==DOUBLE) {
typeName="double";
}
fprintf(outfile,"%s *_%s_main_old = new %s[_%s_main_ncomponents];\n",typeName,fieldName,typeName,fieldName);
fprintf(outfile,"\n");
for(long unsigned int 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");
list<const xmdsVector*> crossVectorList;
for(list<XMLString>::const_iterator pXMLString = crossVectorNamesList()->begin(); pXMLString != crossVectorNamesList()->end(); pXMLString++) {
const xmdsVector* crossVector;
if(!simulation()->field()->getVector(*pXMLString,crossVector)) {
throw xmdsException("Internal error in xmdsIntegrateSI::writeDefines: cannot find cross vector");
}
crossVectorList.push_back(crossVector);
if(crossVector->vectorType()==COMPLEX) {
typeName="complex";
}
else if(crossVector->vectorType()==DOUBLE) {
typeName="double";
}
fprintf(outfile,"%s *_%s_%s_old = new %s[_%s_%s_ncomponents];\n",
typeName,fieldName,crossVector->name()->c_str(),typeName,fieldName,crossVector->name()->c_str());
fprintf(outfile,"\n");
for(long unsigned int i=0;i<crossVector->nComponents();i++) {
fprintf(outfile,"%s d%s_d%s;\n",
typeName,crossVector->componentName(i)->c_str(),
simulation()->field()->geometry()->dimension(crossDimNumber())->name.c_str());
}
fprintf(outfile,"\n");
}
if((simulation()->parameters()->stochastic)&&(!noNoises())) {
fprintf(outfile,"double *_noises = new double[_n_noises];\n");
if(simulation()->parameters()->errorCheck) {
fprintf(outfile,"double *_noises2 = new double[_n_noises];\n");
}
fprintf(outfile,"const double _var = 1/_step");
for(long unsigned int i=0;i<nDims;i++) {
fprintf(outfile,"/_%s_dx%li",fieldName,i);
}
fprintf(outfile,";\n");
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,"");
list<XMLString> myNonLoopPropagationCodeList=*nonLoopPropagationCodeList();
list<integrateMomentGroup> myIntegrateMomentGroupList=*integrateMomentGroupList();
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 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");
// 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,"\n");
char indent[64];
for(long unsigned int i=0;i<nDims;i++) {
indent[i]=0x09;
}
indent[nDims]=0;
if((simulation()->parameters()->stochastic)&&(!noNoises())) {
if(simulation()->parameters()->errorCheck) {
fprintf(outfile,"%sif(_generator_flag==1)\n",indent);
fprintf(outfile,"%s _make_noises(_gen1,_var,_noises,_n_noises);\n",indent);
fprintf(outfile,"%selse if(_generator_flag==2)\n",indent);
fprintf(outfile,"%s _make_noises(_gen2,_var,_noises,_n_noises);\n",indent);
fprintf(outfile,"%selse {\n",indent);
fprintf(outfile,"\n");
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,"\n");
fprintf(outfile,"%s for(unsigned long _s0=0;_s0<_n_noises;_s0++)\n",indent);
fprintf(outfile,"%s _noises[_s0] += _noises2[_s0];\n",indent);
fprintf(outfile,"%s }\n",indent);
fprintf(outfile,"\n");
}
else {
fprintf(outfile,"%s_make_noises(_gen,_var,_noises,_n_noises);\n",indent);
fprintf(outfile,"\n");
}
}
fprintf(outfile,"%sfor(unsigned long _i%li=0;_i%li<_%s_main_ncomponents;_i%li++)\n",indent,nDims+1,nDims+1,fieldName,nDims+1);
fprintf(outfile,"%s _%s_main_old[_i%li] = _active_%s_main[_%s_main_index_pointer + _i%li];\n",indent,fieldName,nDims+1,fieldName,fieldName,nDims+1);
fprintf(outfile,"\n");
if(crossVectorList.size() > 0) {
for(list<const xmdsVector*>::const_iterator ppxmdsVector = crossVectorList.begin(); ppxmdsVector != crossVectorList.end(); ppxmdsVector++) {
const char* cvName = (*ppxmdsVector)->name()->c_str();
fprintf(outfile,"%sfor(unsigned long _i%li=0;_i%li<_%s_%s_ncomponents;_i%li++)\n",
indent,nDims+1,nDims+1,fieldName,cvName,nDims+1);
fprintf(outfile,"%s _%s_%s_old[_i%li] = _active_%s_%s[_%s_%s_index_pointer + _i%li];\n" ,
indent,fieldName,cvName,nDims+1,fieldName,cvName,fieldName,cvName,nDims+1);
fprintf(outfile,"\n");
}
}
if(nIterations()>1) {
fprintf(outfile,"%sfor(unsigned long _i%li=0;_i%li<%li;_i%li++) {\n",indent,nDims+1,nDims+1,nIterations()-1,nDims+1);
fprintf(outfile,"\n");
fprintf(outfile,"// *** propagation (and cross_propagation code) ***\n");
fprintf(outfile,"%s\n",propagationCode()->c_str());
if(crossVectorList.size() > 0) {
fprintf(outfile,"%s\n",crossPropagationCode()->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] = _%s_main_old[%li] + d%s_d%s*(_step/2);\n",
indent,fieldName,fieldName,i,fieldName,i,
mainVector->componentName(i)->c_str(),simulation()->parameters()->propDimName.c_str());
}
for(list<const xmdsVector*>::const_iterator ppxmdsVector = crossVectorList.begin(); ppxmdsVector != crossVectorList.end(); ppxmdsVector++) {
const char* cvName = (*ppxmdsVector)->name()->c_str();
fprintf(outfile,"\n");
for(long unsigned int i=0;i<(*ppxmdsVector)->nComponents();i++) {
fprintf(outfile,"%s _active_%s_%s[_%s_%s_index_pointer + %li] = _%s_%s_old[%li] + d%s_d%s*(_%s_dx%li/2);\n",
indent,fieldName,cvName,fieldName,cvName,i,fieldName,cvName,i,(*ppxmdsVector)->componentName(i)->c_str(),
simulation()->field()->geometry()->dimension(crossDimNumber())->name.c_str(),fieldName,crossDimNumber());
}
}
fprintf(outfile,"%s }\n",indent);
fprintf(outfile,"\n");
}
fprintf(outfile,"// *** propagation (and cross_propagation code) ***\n");
fprintf(outfile,"%s\n",propagationCode()->c_str());
if(crossVectorList.size() > 0) {
fprintf(outfile,"%s\n",crossPropagationCode()->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] = _%s_main_old[%li] + d%s_d%s*_step;\n",
indent,fieldName,fieldName,i,fieldName,i,
mainVector->componentName(i)->c_str(),simulation()->parameters()->propDimName.c_str());
}
if(crossVectorList.size() > 0) {
fprintf(outfile,"\n");
fprintf(outfile,"%sif(_i%li<_main_lattice%li-1) {",indent,crossDimNumber(),crossDimNumber());
}
for(list<const xmdsVector*>::const_iterator ppxmdsVector = crossVectorList.begin(); ppxmdsVector != crossVectorList.end(); ppxmdsVector++) {
const char* cvName = (*ppxmdsVector)->name()->c_str();
fprintf(outfile,"\n");
for(long unsigned int i=0;i<(*ppxmdsVector)->nComponents();i++) {
fprintf(outfile,"%s_active_%s_%s[_%s_%s_index_pointer + %li + _%s_%s_ncomponents",indent,fieldName,cvName,fieldName,cvName,i,fieldName,cvName);
for(long unsigned int i2=crossDimNumber()+1;i2<nDims;i2++) {
fprintf(outfile,"*_%s_lattice%li",fieldName,i2);
}
fprintf(outfile,"] = _%s_%s_old[%li] + d%s_d%s*_%s_dx%li;\n",
fieldName,cvName,i,(*ppxmdsVector)->componentName(i)->c_str(),
simulation()->field()->geometry()->dimension(crossDimNumber())->name.c_str(),fieldName,crossDimNumber());
}
}
if(crossVectorList.size() > 0) {
fprintf(outfile,"%s }\n\n",indent);
for(list<const xmdsVector*>::const_iterator ppxmdsVector = crossVectorList.begin(); ppxmdsVector != crossVectorList.end(); ppxmdsVector++) {
const char* cvName = (*ppxmdsVector)->name()->c_str();
fprintf(outfile,"\n");
fprintf(outfile,"%s for(unsigned long _i%li=0;_i%li<_%s_%s_ncomponents;_i%li++) {\n",
indent,nDims+1,nDims+1,fieldName,cvName,nDims+1);
fprintf(outfile,"%s _active_%s_%s[_%s_%s_index_pointer + _i%li] = _%s_%s_old[_i%li];\n",
indent,fieldName,cvName,fieldName,cvName,nDims+1,fieldName,cvName,nDims+1);
}
fprintf(outfile,"%s }\n",indent);
}
simulation()->field()->closeLoops(outfile,0,myTotalVectorsList);
}
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);
}
// clean up mallocs
fprintf(outfile," delete[] _%s_main_old;\n",fieldName);
for(list<XMLString>::const_iterator pXMLString = crossVectorNamesList()->begin(); pXMLString != crossVectorNamesList()->end(); pXMLString++) {
fprintf(outfile," delete[] _%s_%s_old;\n",
fieldName,pXMLString->c_str());
}
if((simulation()->parameters()->stochastic)&&(!noNoises())) {
fprintf(outfile," delete[] _noises;\n");
if(simulation()->parameters()->errorCheck) {
fprintf(outfile," delete[] _noises2;\n");
}
fprintf(outfile,"\n");
}
fprintf(outfile,"}\n");
fprintf(outfile,"\n");
};
// ******************************************************************************
void xmdsIntegrateSIIP::writeSingleStepCode(
FILE *const outfile,
const stepCaseEnum& stepCase) const {
if(debugFlag) {
printf("xmdsIntegrateSIIP::writeSingleStepCode\n");
}
const char *const propDim = simulation()->parameters()->propDimName.c_str();
const char* kStep="";
if(!constantK()) {
kStep="_step/2";
}
const char* generatorFlag="";
if((simulation()->parameters()->stochastic)&&(!noNoises())) {
if(stepCase==FULLSTEP) {
generatorFlag=",0";
}
else if(stepCase==FIRST_HALFSTEP) {
generatorFlag=",1";
}
else if(stepCase==SECOND_HALFSTEP) {
generatorFlag=",2";
}
}
const char* indent = " ";
if(simulation()->parameters()->errorCheck) {
indent = " ";
}
if(usesKOperators()) {
fprintf(outfile,"%s_segment%li_k_propagate(%s);\n",indent,segmentNumber,kStep);
fprintf(outfile,"\n");
}
fprintf(outfile,"%s%s += _step/2;\n",indent,propDim);
fprintf(outfile,"\n");
fprintf(outfile,"%s_segment%li_x_propagate(_step,cycle%s);\n",indent,segmentNumber,generatorFlag);
fprintf(outfile,"\n");
if(usesKOperators()) {
fprintf(outfile,"%s_segment%li_k_propagate(%s);\n",indent,segmentNumber,kStep);
fprintf(outfile,"\n");
}
fprintf(outfile,"%s%s += _step/2;\n",indent,propDim);
fprintf(outfile,"\n");
};
syntax highlighted by Code2HTML, v. 0.9.1