/*
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: xmdsintegrateip.cc,v 1.14 2005/05/02 02:35:58 sebwuester Exp $
*/
/*! @file xmdsintegrateip.cc
@brief Integrate element parsing classes and methods; interaction picture
More detailed explanation...
*/
#include<xmlbasics.h>
#include<dom3.h>
#include<xmdsutils.h>
#include<xmdsclasses.h>
// ******************************************************************************
// ******************************************************************************
// xmdsIntegrateIP public
// ******************************************************************************
// ******************************************************************************
extern bool debugFlag;
long nxmdsIntegrateIPs=0; //!< Number of integrate IP objects
// ******************************************************************************
xmdsIntegrateIP::xmdsIntegrateIP(
const xmdsSimulation *const yourSimulation,
const bool& yourVerboseMode) :
xmdsIntegrate(yourSimulation,yourVerboseMode) {
if(debugFlag) {
nxmdsIntegrateIPs++;
printf("xmdsIntegrateIP::xmdsIntegrateIP\n");
printf("nxmdsIntegrateIPs=%li\n",nxmdsIntegrateIPs);
}
};
// ******************************************************************************
xmdsIntegrateIP::~xmdsIntegrateIP() {
if(debugFlag) {
nxmdsIntegrateIPs--;
printf("xmdsIntegrateIP::~xmdsIntegrateIP\n");
printf("nxmdsIntegrateIPs=%li\n",nxmdsIntegrateIPs);
}
};
// ******************************************************************************
void xmdsIntegrateIP::processElement(
const Element *const yourElement) {
if(debugFlag) {
printf("xmdsIntegrateIP::processElement\n");
}
// ************************************
// parse code for operators
const xmdsVector* mainVector;
if(!simulation()->field()->getVector("main",mainVector)) {
throw xmdsException("Internal error in xmdsIntegrateIP::processElement: cannot find 'main' vector");
}
XMLString* theCode = propagationCode();
XMLString nextOperatorName;
XMLString nextComponentName;
unsigned long start=0;
unsigned long end=0;
while(end<theCode->length()) {
if(findNextcoPair(nextOperatorName,nextComponentName,start,end)) {
unsigned long nextComponentNumber;
if(!mainVector->getComponent(nextComponentName,nextComponentNumber)) {
sprintf(errorMessage(),"[%s] is not a component of the main vector",nextComponentName.c_str());
throw xmdsException(yourElement,errorMessage());
}
unsigned long nextOperatorNumber;
if(!getKOperator(nextOperatorName,nextOperatorNumber)) {
sprintf(errorMessage(),"'%s' was not defined in <k_operators>",nextOperatorName.c_str());
throw xmdsException(yourElement,errorMessage());
}
unsigned long dummycoKey;
if(getcoKey(nextComponentNumber,nextOperatorNumber,dummycoKey)) {
sprintf(errorMessage(),"'%s[%s]' used multiply",nextOperatorName.c_str(),nextComponentName.c_str());
throw xmdsException(yourElement,errorMessage());
}
if(verbose()) {
printf("adding operator-component pair: %s[%s]\n",nextOperatorName.c_str(),nextComponentName.c_str());
}
addcoPair(nextComponentNumber,nextOperatorNumber);
theCode->replaceData(start,end-start,"0");
start = start+1;
}
else {
start = end;
}
}
};
// ******************************************************************************
// ******************************************************************************
// xmdsIntegrateIP protected
// ******************************************************************************
// *****************************************************************************
// ******************************************************************************
void xmdsIntegrateIP::writePrototypes(
FILE *const outfile) const {
if(debugFlag) {
printf("xmdsIntegrateIP::writePrototypes\n");
}
if(usesKOperators()) {
if(AdaptiveIP()){
if(Smallmemory()){
fprintf(outfile,"void _segment%li_k_propagate(\n",segmentNumber);
fprintf(outfile," const double& _step);\n");
fprintf(outfile,"\n");
}
else{
fprintf(outfile,"void _segment%li_calculate_k_operator_field(double _step);\n",segmentNumber);
fprintf(outfile,"\n");
fprintf(outfile,"void _segment%li_k_propagate(int exponent);\n",segmentNumber);
fprintf(outfile,"\n");
}
}else{
if(constantK()) {
fprintf(outfile,"void _segment%li_calculate_k_operator_field();\n",segmentNumber);
fprintf(outfile,"\n");
fprintf(outfile,"void _segment%li_k_propagate();\n",segmentNumber);
fprintf(outfile,"\n");
}
else {
fprintf(outfile,"void _segment%li_k_propagate(\n",segmentNumber);
fprintf(outfile," const double& _step);\n");
fprintf(outfile,"\n");
}
}
}
};
// ******************************************************************************
void xmdsIntegrateIP::writeRoutines(
FILE *const outfile) const {
if(debugFlag) {
printf("xmdsIntegrateIP::writeRoutines\n");
}
if(usesKOperators()) {
if(AdaptiveIP()){
if(Smallmemory())
writeTimeDepkPropagateRoutine(outfile); // use conventional routine
else{
if(constantK()) // optimized for speed, but needs lots of memory
writeARK45CalculatekOperatorFieldRoutine(outfile);
else // optimized for speed, but needs even more memory
writeARK45TimeDepkCalculatekOperatorFieldRoutine(outfile);
writeARK45kPropagateRoutine(outfile);
}
}else{
if(constantK()) {
writeCalculatekOperatorFieldRoutine(outfile);
writeConstkPropagateRoutine(outfile);
}
else
writeTimeDepkPropagateRoutine(outfile);
}
}
};
// ******************************************************************************
// ******************************************************************************
// xmdsIntegrateIP private
// ******************************************************************************
// ******************************************************************************
// ******************************************************************************
void xmdsIntegrateIP::writeCalculatekOperatorFieldRoutine(
FILE *const outfile) const {
if(debugFlag) {
printf("xmdsIntegrateIP::writeCalculatekOperatorFieldRoutine\n");
}
const unsigned long nDims = simulation()->field()->geometry()->nDims();
const unsigned long fullSpace = simulation()->field()->geometry()->fullSpace();
unsigned long i;
fprintf(outfile,"// *************************\n");
fprintf(outfile,"void _segment%li_calculate_k_operator_field() {\n",segmentNumber);
fprintf(outfile,"\n");
for(i=0;i<nKOperators();i++) {
fprintf(outfile,"complex %s;\n",KOperator(i)->c_str());
}
fprintf(outfile,"\n");
if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
fprintf(outfile,"_segment%li_k_operator_field = new complex[total_local_size*_segment%li_nkoperators];\n",
segmentNumber,segmentNumber);
fprintf(outfile,"\n");
}
fprintf(outfile,"double _step=(%s/(double)%li)/2;\n",interval()->c_str(),lattice());
fprintf(outfile,"\n");
if(simulation()->parameters()->errorCheck) {
fprintf(outfile,"if(_half_step)\n");
fprintf(outfile," _step=_step/2;\n");
fprintf(outfile,"\n");
}
simulation()->field()->vectors2space(outfile,fullSpace,*KVectorNamesList(),"");
fprintf(outfile,"unsigned long _k_operator_index_pointer=0;\n");
simulation()->field()->openLoops(outfile,fullSpace,*KVectorNamesList());
char indent[64];
for(i=0;i<nDims;i++) {
indent[i]=0x09;
}
indent[nDims]=0;
fprintf(outfile,"\n");
fprintf(outfile,"// ********** Code from k_operators *************\n");
fprintf(outfile,"%s\n",KOperatorsCode()->c_str());
fprintf(outfile,"// **********************************************\n");
fprintf(outfile,"\n");
for(i=0;i<nKOperators();i++) {
fprintf(outfile,"%s_segment%li_k_operator_field[_k_operator_index_pointer + %li].re = exp(%s.re*_step)*cos(%s.im*_step);\n",
indent,segmentNumber,i,KOperator(i)->c_str(),KOperator(i)->c_str());
fprintf(outfile,"%s_segment%li_k_operator_field[_k_operator_index_pointer + %li].im = exp(%s.re*_step)*sin(%s.im*_step);\n",
indent,segmentNumber,i,KOperator(i)->c_str(),KOperator(i)->c_str());
}
fprintf(outfile,"\n");
fprintf(outfile,"%s_k_operator_index_pointer += _segment%li_nkoperators;\n",indent,segmentNumber);
simulation()->field()->closeLoops(outfile,fullSpace,*KVectorNamesList());
fprintf(outfile,"}\n");
fprintf(outfile,"\n");
};
// ******************************************************************************
void xmdsIntegrateIP::writeConstkPropagateRoutine(
FILE *const outfile) const {
if(debugFlag) {
printf("xmdsIntegrateIP::writeConstkPropagateRoutine\n");
}
const unsigned long fullSpace = simulation()->field()->geometry()->fullSpace();
const char *const fieldName = simulation()->field()->name()->c_str();
list<XMLString> tempVectorNamesList;
tempVectorNamesList.push_back("main");
const xmdsVector* mainVector;
if(!simulation()->field()->getVector("main",mainVector)) {
throw xmdsException("Internal error in xmdsIntegrateIP::writeConstkPropagateRoutine: cannot find 'main' vector");
}
fprintf(outfile,"// *************************\n");
fprintf(outfile,"void _segment%li_k_propagate() {\n",segmentNumber);
fprintf(outfile,"\n");
fprintf(outfile,"double _temp;\n");
fprintf(outfile,"unsigned long _segment%li_kop_index_pointer=0;\n",segmentNumber);
fprintf(outfile,"unsigned long _active_%s_index_pointer=0;\n",fieldName);
fprintf(outfile,"\n");
simulation()->field()->vectors2space(outfile,fullSpace,tempVectorNamesList,"");
if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
fprintf(outfile,"for(long _i0=0;_i0<total_local_size;_i0++) {\n");
}
else {
fprintf(outfile,"for(unsigned long _i0=0;_i0<_%s_size;_i0++) {\n",fieldName);
}
fprintf(outfile,"\n");
for(unsigned long i=0; i<mainVector->nComponents(); i++) {
const coStruct* thecoStruct;
if(getcoStruct(i,thecoStruct)) {
for(list<unsigned long>::const_iterator pULong = thecoStruct->operatorNumbersList.begin(); pULong != thecoStruct->operatorNumbersList.end(); pULong++) {
fprintf(outfile," _temp = _segment%li_k_operator_field[_segment%li_kop_index_pointer + %li].re\n",
segmentNumber,segmentNumber,*pULong);
fprintf(outfile," *_active_%s_main[_active_%s_index_pointer + %li].re\n",fieldName,fieldName,i);
fprintf(outfile," - _segment%li_k_operator_field[_segment%li_kop_index_pointer + %li].im\n",
segmentNumber,segmentNumber,*pULong);
fprintf(outfile," *_active_%s_main[_active_%s_index_pointer + %li].im;\n",fieldName,fieldName,i);
fprintf(outfile,"\n");
fprintf(outfile," _active_%s_main[_active_%s_index_pointer + %li].im =\n",fieldName,fieldName,i);
fprintf(outfile," _segment%li_k_operator_field[_segment%li_kop_index_pointer + %li].re\n",
segmentNumber,segmentNumber,*pULong);
fprintf(outfile," *_active_%s_main[_active_%s_index_pointer + %li].im\n",fieldName,fieldName,i);
fprintf(outfile," + _segment%li_k_operator_field[_segment%li_kop_index_pointer + %li].im\n",
segmentNumber,segmentNumber,*pULong);
fprintf(outfile," *_active_%s_main[_active_%s_index_pointer + %li].re;\n",fieldName,fieldName,i);
fprintf(outfile,"\n");
fprintf(outfile," _active_%s_main[_active_%s_index_pointer + %li].re=_temp;\n",fieldName,fieldName,i);
fprintf(outfile,"\n");
}
}
}
fprintf(outfile," _segment%li_kop_index_pointer+=_segment%li_nkoperators;\n",segmentNumber,segmentNumber);
fprintf(outfile," _active_%s_index_pointer+=_%s_main_ncomponents;\n",fieldName,fieldName);
fprintf(outfile," };\n");
fprintf(outfile,"}\n");
fprintf(outfile,"\n");
};
// ******************************************************************************
void xmdsIntegrateIP::writeTimeDepkPropagateRoutine(
FILE *const outfile) const {
if(debugFlag) {
printf("xmdsIntegrateIP::writeTimeDepkPropagateRoutine\n");
}
const unsigned long nDims = simulation()->field()->geometry()->nDims();
const unsigned long fullSpace = simulation()->field()->geometry()->fullSpace();
const char *const fieldName = simulation()->field()->name()->c_str();
const xmdsVector* mainVector;
if(!simulation()->field()->getVector("main",mainVector)) {
throw xmdsException("Internal error in xmdsIntegrateIP::writeTimeDepkPropagateRoutine: cannot find 'main' vector");
}
fprintf(outfile,"// *************************\n");
fprintf(outfile,"void _segment%li_k_propagate(\n",segmentNumber);
fprintf(outfile," const double& _step) {\n");
fprintf(outfile,"\n");
for(unsigned long i=0;i<nKOperators();i++) {
fprintf(outfile,"complex %s;\n",KOperator(i)->c_str());
}
fprintf(outfile,"\n");
fprintf(outfile,"double _temp;\n");
fprintf(outfile,"double _temp2;\n");
fprintf(outfile,"\n");
simulation()->field()->vectors2space(outfile,fullSpace,*KVectorNamesList(),"");
simulation()->field()->openLoops(outfile,fullSpace,*KVectorNamesList());
char indent[64];
for(unsigned long i=0;i<nDims;i++) {
indent[i]=0x09;
}
indent[nDims]=0;
fprintf(outfile,"\n");
fprintf(outfile,"// ********** Code from k_operators *************\n");
fprintf(outfile,"%s\n",KOperatorsCode()->c_str());
fprintf(outfile,"// **********************************************\n");
fprintf(outfile,"\n");
for(unsigned long i=0;i<nKOperators();i++) {
fprintf(outfile,"%s_temp2 = exp(%s.re*_step);\n",indent,KOperator(i)->c_str());
fprintf(outfile,"%s_temp = _temp2*cos(%s.im*_step);\n",indent,KOperator(i)->c_str());
fprintf(outfile,"%s%s.im = _temp2*sin(%s.im*_step);\n",indent,KOperator(i)->c_str(),KOperator(i)->c_str());
fprintf(outfile,"%s%s.re = _temp;\n",indent,KOperator(i)->c_str());
fprintf(outfile,"\n");
}
for(unsigned long i=0; i<mainVector->nComponents(); i++) {
const coStruct* thecoStruct;
if(getcoStruct(i,thecoStruct)) {
for(list<unsigned long>::const_iterator pULong = thecoStruct->operatorNumbersList.begin(); pULong != thecoStruct->operatorNumbersList.end(); pULong++) {
fprintf(outfile,"%s_temp = %s.re*_active_%s_main[_%s_main_index_pointer + %li].re\n",indent,KOperator(*pULong)->c_str(),fieldName,fieldName,i);
fprintf(outfile,"%s - %s.im*_active_%s_main[_%s_main_index_pointer + %li].im;\n",indent,KOperator(*pULong)->c_str(),fieldName,fieldName,i);
fprintf(outfile,"\n");
fprintf(outfile,"%s_active_%s_main[_%s_main_index_pointer + %li].im = %s.re*_active_%s_main[_%s_main_index_pointer + %li].im\n",
indent,fieldName,fieldName,i,KOperator(*pULong)->c_str(),fieldName,fieldName,i);
fprintf(outfile,"%s + %s.im*_active_%s_main[_%s_main_index_pointer + %li].re;\n",
indent,KOperator(*pULong)->c_str(),fieldName,fieldName,i);
fprintf(outfile,"\n");
fprintf(outfile,"%s_active_%s_main[_%s_main_index_pointer + %li].re=_temp;\n",indent,fieldName,fieldName,i);
fprintf(outfile,"\n");
}
}
}
simulation()->field()->closeLoops(outfile,fullSpace,*KVectorNamesList());
fprintf(outfile,"}\n");
fprintf(outfile,"\n");
};
// ******************************************************************************
void xmdsIntegrateIP::writeARK45CalculatekOperatorFieldRoutine(
FILE *const outfile) const {
if(debugFlag) {
printf("xmdsIntegrateARK45IP::writeCalculatekOperatorFieldRoutine\n");
}
const unsigned long nDims = simulation()->field()->geometry()->nDims();
const unsigned long fullSpace = simulation()->field()->geometry()->fullSpace();
unsigned long i;
fprintf(outfile,"// *************************\n");
fprintf(outfile,"void _segment%li_calculate_k_operator_field(double _step) {\n",segmentNumber);
fprintf(outfile,"\n");
for(i=0;i<nKOperators();i++) {
fprintf(outfile,"complex %s;\n",KOperator(i)->c_str());
}
fprintf(outfile,"\n");
fprintf(outfile,"double factor1=1.0/5; // a_raw[2]\n");
fprintf(outfile,"double factor2=1.0/10; // a_raw[3]-a_raw[2]\n");
fprintf(outfile,"double factor3=2.0/5; // a_raw[4]-a_raw[2]\n");
fprintf(outfile,"double factor4=4.0/5; // a_raw[5]-a_raw[2]\n");
fprintf(outfile,"double factor5=27.0/40; // a_raw[6]-a_raw[2]\n");
simulation()->field()->vectors2space(outfile,fullSpace,*KVectorNamesList(),"");
fprintf(outfile,"unsigned long _k_operator_index_pointer=0;\n");
simulation()->field()->openLoops(outfile,fullSpace,*KVectorNamesList());
char indent[64];
for(i=0;i<nDims;i++) {
indent[i]=0x09;
}
indent[nDims]=0;
fprintf(outfile,"\n");
fprintf(outfile,"// ********** Code from k_operators *************\n");
fprintf(outfile,"%s\n",KOperatorsCode()->c_str());
fprintf(outfile,"// **********************************************\n");
fprintf(outfile,"\n");
for(int n=1;n<6;n++){
for(i=0;i<nKOperators();i++) {
fprintf(outfile,"%s_segment%li_k_operator_field_%i[_k_operator_index_pointer + %li].re = exp(%s.re*factor%i*_step)*cos(%s.im*factor%i*_step);\n",
indent,segmentNumber,n,i,KOperator(i)->c_str(),n,KOperator(i)->c_str(),n);
fprintf(outfile,"%s_segment%li_k_operator_field_%i[_k_operator_index_pointer + %li].im = exp(%s.re*factor%i*_step)*sin(%s.im*factor%i*_step);\n",
indent,segmentNumber,n,i,KOperator(i)->c_str(),n,KOperator(i)->c_str(),n);
}
fprintf(outfile,"\n");
}
fprintf(outfile,"\n");
fprintf(outfile,"%s_k_operator_index_pointer += _segment%li_nkoperators;\n",indent,segmentNumber);
simulation()->field()->closeLoops(outfile,fullSpace,*KVectorNamesList());
fprintf(outfile,"}\n");
fprintf(outfile,"\n");
};
// ******************************************************************************
void xmdsIntegrateIP::writeARK45kPropagateRoutine(
FILE *const outfile) const {
if(debugFlag) {
printf("xmdsIntegrateARK45IP::writeConstkPropagateRoutine\n");
}
const unsigned long fullSpace = simulation()->field()->geometry()->fullSpace();
const char *const fieldName = simulation()->field()->name()->c_str();
list<XMLString> tempVectorNamesList;
tempVectorNamesList.push_back("main");
const xmdsVector* mainVector;
if(!simulation()->field()->getVector("main",mainVector)) {
throw xmdsException("Internal error in xmdsIntegrateARK45IP::writeConstkPropagateRoutine: cannot find 'main' vector");
}
fprintf(outfile,"// *************************\n");
fprintf(outfile,"void _segment%li_k_propagate(int exponent) {\n",segmentNumber);
fprintf(outfile,"\n");
fprintf(outfile,"double _temp;\n");
fprintf(outfile,"unsigned long _segment%li_kop_index_pointer=0;\n",segmentNumber);
fprintf(outfile,"unsigned long _active_%s_index_pointer=0;\n",fieldName);
fprintf(outfile,"\n");
simulation()->field()->vectors2space(outfile,fullSpace,tempVectorNamesList,"");
fprintf(outfile,"if(exponent>0){\n");
fprintf(outfile,"switch(exponent){\n");
for(int n=1;n<6;n++){
fprintf(outfile,"case %i:\n",n);
if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
fprintf(outfile,"for(long _i0=0;_i0<total_local_size;_i0++) {\n");
}
else {
fprintf(outfile,"for(unsigned long _i0=0;_i0<_%s_size;_i0++) {\n",fieldName);
}
fprintf(outfile,"\n");
for(unsigned long i=0; i<mainVector->nComponents(); i++) {
const coStruct* thecoStruct;
if(getcoStruct(i,thecoStruct)) {
for(list<unsigned long>::const_iterator pULong = thecoStruct->operatorNumbersList.begin(); pULong != thecoStruct->operatorNumbersList.end(); pULong++) {
fprintf(outfile," _temp = _segment%li_k_operator_field_%i[_segment%li_kop_index_pointer + %li].re\n",
segmentNumber,n,segmentNumber,*pULong);
fprintf(outfile," *_active_%s_main[_active_%s_index_pointer + %li].re\n",fieldName,fieldName,i);
fprintf(outfile," - _segment%li_k_operator_field_%i[_segment%li_kop_index_pointer + %li].im\n",
segmentNumber,n,segmentNumber,*pULong);
fprintf(outfile," *_active_%s_main[_active_%s_index_pointer + %li].im;\n",fieldName,fieldName,i);
fprintf(outfile,"\n");
fprintf(outfile," _active_%s_main[_active_%s_index_pointer + %li].im =\n",fieldName,fieldName,i);
fprintf(outfile," _segment%li_k_operator_field_%i[_segment%li_kop_index_pointer + %li].re\n",
segmentNumber,n,segmentNumber,*pULong);
fprintf(outfile," *_active_%s_main[_active_%s_index_pointer + %li].im\n",fieldName,fieldName,i);
fprintf(outfile," + _segment%li_k_operator_field_%i[_segment%li_kop_index_pointer + %li].im\n",
segmentNumber,n,segmentNumber,*pULong);
fprintf(outfile," *_active_%s_main[_active_%s_index_pointer + %li].re;\n",fieldName,fieldName,i);
fprintf(outfile,"\n");
fprintf(outfile," _active_%s_main[_active_%s_index_pointer + %li].re=_temp;\n",fieldName,fieldName,i);
fprintf(outfile,"\n");
}
}
}
fprintf(outfile," _segment%li_kop_index_pointer+=_segment%li_nkoperators;\n",segmentNumber,segmentNumber);
fprintf(outfile," _active_%s_index_pointer+=_%s_main_ncomponents;\n",fieldName,fieldName);
fprintf(outfile," };\n");
fprintf(outfile," break;\n");
}
fprintf(outfile,"}\n");
fprintf(outfile,"} else {\n");
fprintf(outfile,"exponent=abs(exponent);\n");
fprintf(outfile,"switch(exponent){\n");
for(int n=1;n<6;n++){
fprintf(outfile,"case %i:\n",n);
if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
fprintf(outfile,"for(long _i0=0;_i0<total_local_size;_i0++) {\n");
}
else {
fprintf(outfile,"for(unsigned long _i0=0;_i0<_%s_size;_i0++) {\n",fieldName);
}
fprintf(outfile,"\n");
for(unsigned long i=0; i<mainVector->nComponents(); i++) {
const coStruct* thecoStruct;
if(getcoStruct(i,thecoStruct)) {
for(list<unsigned long>::const_iterator pULong = thecoStruct->operatorNumbersList.begin(); pULong != thecoStruct->operatorNumbersList.end(); pULong++) {
fprintf(outfile," _active_%s_main[_active_%s_index_pointer + %li]/=_segment%li_k_operator_field_%i[_segment%li_kop_index_pointer + %li];\n",fieldName,fieldName,i,segmentNumber,n,segmentNumber,*pULong);
fprintf(outfile,"\n");
}
}
}
fprintf(outfile," _segment%li_kop_index_pointer+=_segment%li_nkoperators;\n",segmentNumber,segmentNumber);
fprintf(outfile," _active_%s_index_pointer+=_%s_main_ncomponents;\n",fieldName,fieldName);
fprintf(outfile," };\n");
fprintf(outfile," break;\n");
}
fprintf(outfile,"}\n");
fprintf(outfile,"}\n");
fprintf(outfile,"}\n");
fprintf(outfile,"\n");
};
// ******************************************************************************
void xmdsIntegrateIP::writeARK45TimeDepkCalculatekOperatorFieldRoutine(
FILE *const outfile) const {
if(debugFlag) {
printf("xmdsIntegrateARK45IP::writeARK45TimeDepkCalculatekOperatorFieldRoutine\n");
}
const unsigned long nDims = simulation()->field()->geometry()->nDims();
const char *const propDim = simulation()->parameters()->propDimName.c_str();
const unsigned long fullSpace = simulation()->field()->geometry()->fullSpace();
unsigned long i;
fprintf(outfile,"// *************************\n");
fprintf(outfile,"void _segment%li_calculate_k_operator_field(double _step) {\n",segmentNumber);
fprintf(outfile,"\n");
for(i=0;i<nKOperators();i++) {
fprintf(outfile,"complex %s;\n",KOperator(i)->c_str());
}
fprintf(outfile,"\n");
fprintf(outfile,"double factor1=1.0/5; // a_raw[2]\n");
fprintf(outfile,"double factor2=1.0/10; // a_raw[3]-a_raw[2]\n");
fprintf(outfile,"double factor3=2.0/5; // a_raw[4]-a_raw[2]\n");
fprintf(outfile,"double factor4=4.0/5; // a_raw[5]-a_raw[2]\n");
fprintf(outfile,"double factor5=27.0/40; // a_raw[6]-a_raw[2]\n");
fprintf(outfile,"double _store_time=t;\n");
fprintf(outfile,"\n");
fprintf(outfile,"double _store_%s=%s; \n",propDim,propDim);
fprintf(outfile,"double _d%s1=0.0; \n",propDim);
fprintf(outfile,"double _d%s2=3.0/10*_step; \n",propDim);
fprintf(outfile,"double _d%s3=3.0/10*_step; \n",propDim);
fprintf(outfile,"double _d%s4=2.0/5*_step; \n",propDim);
fprintf(outfile,"double _d%s5=-1.0/8*_step; \n",propDim);
fprintf(outfile,"\n");
simulation()->field()->vectors2space(outfile,fullSpace,*KVectorNamesList(),"");
fprintf(outfile,"unsigned long _k_operator_index_pointer;\n");
for(int n=1;n<6;n++){
fprintf(outfile,"_k_operator_index_pointer=0;\n");
fprintf(outfile,"%s+=_d%s%i;\n",propDim,propDim,n);
fprintf(outfile,"\n");
fprintf(outfile,"{\n"); //make counters in loops local
simulation()->field()->openLoops(outfile,fullSpace,*KVectorNamesList());
char indent[64];
for(i=0;i<nDims;i++) {
indent[i]=0x09;
}
indent[nDims]=0;
fprintf(outfile,"\n");
fprintf(outfile,"// ********** Code from k_operators *************\n");
fprintf(outfile,"%s\n",KOperatorsCode()->c_str());
fprintf(outfile,"// **********************************************\n");
fprintf(outfile,"\n");
for(i=0;i<nKOperators();i++) {
fprintf(outfile,"%s_segment%li_k_operator_field_%i[_k_operator_index_pointer + %li].re = exp(%s.re*factor%i*_step)*cos(%s.im*factor%i*_step);\n",
indent,segmentNumber,n,i,KOperator(i)->c_str(),n,KOperator(i)->c_str(),n);
fprintf(outfile,"%s_segment%li_k_operator_field_%i[_k_operator_index_pointer + %li].im = exp(%s.re*factor%i*_step)*sin(%s.im*factor%i*_step);\n",
indent,segmentNumber,n,i,KOperator(i)->c_str(),n,KOperator(i)->c_str(),n);
}
fprintf(outfile,"\n");
fprintf(outfile,"%s_k_operator_index_pointer += _segment%li_nkoperators;\n",indent,segmentNumber);
simulation()->field()->closeLoops(outfile,fullSpace,*KVectorNamesList());
fprintf(outfile,"}\n");// end-make counters in loops local
fprintf(outfile,"\n");
}
fprintf(outfile,"%s=_store_%s;\n",propDim,propDim);
fprintf(outfile,"}\n");
fprintf(outfile,"\n");
};
syntax highlighted by Code2HTML, v. 0.9.1