/*
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: xmdsintegraterk4.cc,v 1.19 2004/10/06 07:52:38 joehope Exp $
*/
/*! @file xmdsintegraterk4.cc
@brief Integrate element parsing classes and methods; fourth order Runge-Kutta
More detailed explanation...
*/
#include<xmlbasics.h>
#include<dom3.h>
#include<xmdsutils.h>
#include<xmdsclasses.h>
// ******************************************************************************
// ******************************************************************************
// xmdsIntegrateRK4 public
// ******************************************************************************
// ******************************************************************************
extern bool debugFlag;
long nxmdsIntegrateRK4s=0; //!< Number of xmds integrate RK4 objects
// ******************************************************************************
xmdsIntegrateRK4::xmdsIntegrateRK4(
const xmdsSimulation *const yourSimulation,
const bool& yourVerboseMode) :
xmdsIntegrate(yourSimulation,yourVerboseMode) {
if(debugFlag) {
nxmdsIntegrateRK4s++;
printf("xmdsIntegrateRK4::xmdsIntegrateRK4\n");
printf("nxmdsIntegrateRK4s=%li\n",nxmdsIntegrateRK4s);
}
};
// ******************************************************************************
xmdsIntegrateRK4::~xmdsIntegrateRK4() {
if(debugFlag) {
nxmdsIntegrateRK4s--;
printf("xmdsIntegrateRK4::~xmdsIntegrateRK4\n");
printf("nxmdsIntegrateRK4s=%li\n",nxmdsIntegrateRK4s);
}
};
// ******************************************************************************
void xmdsIntegrateRK4::processElement(
const Element *const yourElement) {
if(debugFlag) {
printf("xmdsIntegrateRK4::processElement\n");
}
if((simulation()->parameters()->stochastic)&&(!noNoises())) {
printf("\n");
printf("WARNING: RK4 methods may not always yield correct stochastic convergence.\n");
printf("\n");
}
};
// ******************************************************************************
// ******************************************************************************
// xmdsIntegrateRK4 protected
// ******************************************************************************
// ******************************************************************************
// ******************************************************************************
void xmdsIntegrateRK4::writePrototypes(
FILE *const outfile) const {
if(debugFlag) {
printf("xmdsIntegrateRK4::writePrototypes\n");
}
fprintf(outfile,"void _segment%li(unsigned long cycle);\n",segmentNumber);
fprintf(outfile,"\n");
if(crossVectorNamesList()->size() > 0) {
fprintf(outfile,"void _segment%li_calculate_cross_field();\n",segmentNumber);
fprintf(outfile,"\n");
}
};
// ******************************************************************************
void xmdsIntegrateRK4::writeRoutines(
FILE *const outfile) const {
if(debugFlag) {
printf("xmdsIntegrateRK4::writeRoutines\n");
}
writeMainIntegrateRoutine(outfile);
if(crossVectorNamesList()->size() > 0) {
writeCalculateCrossFieldRoutine(outfile);
}
};
// ******************************************************************************
// ******************************************************************************
// xmdsIntegrateRK4 private
// ******************************************************************************
// ******************************************************************************
// ******************************************************************************
void xmdsIntegrateRK4::writeMainIntegrateRoutine(
FILE *const outfile) const {
if(debugFlag) {
printf("xmdsIntegrateRK4::writeMainIntegrateRoutine\n");
}
const char *const fieldName = simulation()->field()->name()->c_str();
const xmdsVector* mainVector;
if(!simulation()->field()->getVector("main",mainVector)) {
throw xmdsException("Internal error in xmdsIntegrateRK4::writeMainIntegrateRoutine: cannot find 'main' vector");
}
const char* typeName="";
if(mainVector->vectorType()==COMPLEX) {
typeName="complex";
}
else if(mainVector->vectorType()==DOUBLE) {
typeName="double";
}
fprintf(outfile,"/* ******************************************** */\n");
fprintf(outfile,"void _segment%li(unsigned long cycle) {\n",segmentNumber);
fprintf(outfile,"\n");
if((simulation()->parameters()->usempi)&!(simulation()->parameters()->stochastic)){
fprintf(outfile,"%s *akfield_%s_main = new %s[total_local_size*_%s_main_ncomponents];\n",typeName,fieldName,typeName,fieldName);
fprintf(outfile,"%s *aifield_%s_main = new %s[total_local_size*_%s_main_ncomponents];\n",typeName,fieldName,typeName,fieldName);
}
else {
fprintf(outfile,"%s *akfield_%s_main = new %s[_%s_size*_%s_main_ncomponents];\n",typeName,fieldName,typeName,fieldName,fieldName);
fprintf(outfile,"%s *aifield_%s_main = new %s[_%s_size*_%s_main_ncomponents];\n",typeName,fieldName,typeName,fieldName,fieldName);
}
fprintf(outfile,"\n");
if((simulation()->parameters()->stochastic)&&(!noNoises())) {
fprintf(outfile,"const double _var = 1");
for(unsigned long i=0;i<simulation()->field()->geometry()->nDims();i++) {
fprintf(outfile,"/_%s_dx%li",fieldName,i);
}
fprintf(outfile,";\n");
fprintf(outfile,"double *_noise_vector = new double[_%s_size*_n_noises];\n",fieldName);
if(simulation()->parameters()->errorCheck) {
fprintf(outfile,"double *_noise_vector2 = new double[_%s_size*_n_noises];\n",fieldName);
}
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");
if((simulation()->parameters()->stochastic)&&(!noNoises())) {
fprintf(outfile," _make_noises(_gen1,_var/_step,_noise_vector,_%s_size*_n_noises);\n",fieldName);
fprintf(outfile,"\n");
}
writeSingleStepCode(outfile,FIRST_HALFSTEP);
fprintf(outfile,"\n");
if((simulation()->parameters()->stochastic)&&(!noNoises())) {
fprintf(outfile," _make_noises(_gen2,_var/_step,_noise_vector,_%s_size*_n_noises);\n",fieldName);
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");
if((simulation()->parameters()->stochastic)&&(!noNoises())) {
fprintf(outfile," _make_noises(_gen1,_var/_step/2,_noise_vector,_%s_size*_n_noises);\n",fieldName);
fprintf(outfile," _make_noises(_gen2,_var/_step/2,_noise_vector2,_%s_size*_n_noises);\n",fieldName);
fprintf(outfile," for(unsigned long _i1=0;_i1<_%s_size*_n_noises;_i1++)\n",fieldName);
fprintf(outfile," _noise_vector[_i1] += _noise_vector2[_i1];\n");
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");
if((simulation()->parameters()->stochastic)&&(!noNoises())) {
fprintf(outfile," _make_noises(_gen,_var/_step,_noise_vector,_%s_size*_n_noises);\n",fieldName);
fprintf(outfile,"\n");
}
writeSingleStepCode(outfile,FULLSTEP);
}
for(unsigned long 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");
if ((simulation()->parameters()->stochastic)&&(!noNoises())) {
fprintf(outfile, " delete[] _noise_vector;\n");
if (simulation()->parameters()->errorCheck) {
fprintf(outfile, " delete[] _noise_vector2;\n");
}
}
fprintf(outfile, " delete[] akfield_%s_main;\n",fieldName);
fprintf(outfile, " delete[] aifield_%s_main;\n",fieldName);
fprintf(outfile,"}\n");
fprintf(outfile,"\n");
};
// ******************************************************************************
void xmdsIntegrateRK4::writeCalculateCrossFieldRoutine(
FILE *const outfile) const {
if(debugFlag) {
printf("xmdsIntegrateRK4IP::writeCalculateCrossFieldRoutine\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_cross_field() {\n",segmentNumber);
fprintf(outfile,"\n");
if(crossDimNumber()+1<nDims) {
fprintf(outfile,"const unsigned long _%s_cross_size = _%s_lattice%li",fieldName,fieldName,crossDimNumber()+1);
for(unsigned long i=crossDimNumber()+2; i<nDims; i++) {
fprintf(outfile,"*_%s_lattice%li",fieldName,i);
}
fprintf(outfile,";\n");
}
else
fprintf(outfile,"const unsigned long _%s_cross_size = 1;\n",fieldName);
fprintf(outfile,"\n");
// need to create a vector list for ONLY main vectors
list<XMLString> myMainVectorNamesList;
for(list<XMLString>::const_iterator pXMLString = vectorNamesList()->begin(); pXMLString != vectorNamesList()->end(); pXMLString++) {
list<XMLString>::const_iterator pXMLString2 = crossVectorNamesList()->begin();
while((pXMLString2 != crossVectorNamesList()->end()) && (*pXMLString2 != *pXMLString)) {
pXMLString2++;
}
if(*pXMLString2 != *pXMLString) {
myMainVectorNamesList.push_back(*pXMLString);
}
}
const char* typeName;
list<const xmdsVector*> mainVectorList;
for(list<XMLString>::const_iterator pXMLString = myMainVectorNamesList.begin(); pXMLString != myMainVectorNamesList.end(); pXMLString++) {
const xmdsVector* mainVector;
if(!simulation()->field()->getVector(*pXMLString,mainVector)) {
throw xmdsException("Internal error in xmdsIntegrateRK4::writeCalculateCrossFieldRoutine: cannot find main vector");
}
mainVectorList.push_back(mainVector);
if(mainVector->vectorType()==DOUBLE) {
typeName="double";
}
else {
typeName="complex";
}
fprintf(outfile,"%s *_%s_%s_old = new %s[_%s_cross_size*_%s_%s_ncomponents];\n",
typeName,fieldName,pXMLString->c_str(),typeName,fieldName,fieldName,pXMLString->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 xmdsIntegrateRK4::writeCalculateCrossFieldRoutine: cannot find cross vector");
}
crossVectorList.push_back(crossVector);
if(crossVector->vectorType()==DOUBLE) {
typeName="double";
}
else {
typeName="complex";
}
fprintf(outfile,"%s *_%s_%s_K = new %s[_%s_cross_size*_%s_%s_ncomponents];\n",
typeName,fieldName,crossVector->name()->c_str(),typeName,fieldName,fieldName,crossVector->name()->c_str());
fprintf(outfile,"%s *_%s_%s_I = new %s[_%s_cross_size*_%s_%s_ncomponents];\n",
typeName,fieldName,crossVector->name()->c_str(),typeName,fieldName,fieldName,crossVector->name()->c_str());
fprintf(outfile,"%s *_%s_%s_d = new %s[_%s_cross_size*_%s_%s_ncomponents];\n",
typeName,fieldName,crossVector->name()->c_str(),typeName,fieldName,fieldName,crossVector->name()->c_str());
fprintf(outfile,"\n");
for(unsigned long 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");
}
// add cross vectors to total vectors to use
list<XMLString> myTotalVectorsList = myMainVectorNamesList;
for(list<XMLString>::const_iterator pXMLString = crossVectorNamesList()->begin(); pXMLString != crossVectorNamesList()->end(); pXMLString++) {
myTotalVectorsList.push_back(*pXMLString);
}
simulation()->field()->vectors2space(outfile,0,myTotalVectorsList,"");
// open outer loops
for(unsigned long i=0; i<crossDimNumber(); i++) {
for(unsigned long j=0; j<i; j++) {
fprintf(outfile," ");
}
fprintf(outfile,"double %s = _%s_xmin%li;\n",simulation()->field()->geometry()->dimension(i)->name.c_str(),fieldName,i);
fprintf(outfile,"\n");
for(unsigned long j=0; j<i; j++) {
fprintf(outfile," ");
}
fprintf(outfile,"for(unsigned long _i%li=0; _i%li<_%s_lattice%li; _i%li++) {\n",i,i,fieldName,i,i);
fprintf(outfile,"\n");
}
for(unsigned long j=0; j<crossDimNumber(); j++) {
fprintf(outfile," ");
}
fprintf(outfile,"double %s = _%s_xmin%li;\n",
simulation()->field()->geometry()->dimension(crossDimNumber())->name.c_str(),fieldName,crossDimNumber());
fprintf(outfile,"\n");
for(unsigned long j=0; j<crossDimNumber(); j++) {
fprintf(outfile," ");
}
fprintf(outfile,"for(unsigned long _i%li=0; _i%li<_%s_lattice%li-1; _i%li++) {\n",
crossDimNumber(),crossDimNumber(),fieldName,crossDimNumber(),crossDimNumber());
fprintf(outfile,"\n");
char indent[64];
for(unsigned long i=0;i<crossDimNumber()+1;i++) {
indent[i]=0x09;
}
indent[crossDimNumber()+1]=0;
char indent2[64];
for(unsigned long i=0;i<nDims;i++) {
indent2[i]=0x09;
}
indent2[nDims]=0;
for(list<XMLString>::const_iterator pXMLString = myTotalVectorsList.begin(); pXMLString != myTotalVectorsList.end(); pXMLString++) {
fprintf(outfile,"%sunsigned long _%s_%s_index_pointer_begin=0;\n",indent,fieldName,pXMLString->c_str());
for(unsigned long i=0; i<crossDimNumber()+1; i++) {
fprintf(outfile,"%s_%s_%s_index_pointer_begin += _i%li",indent,fieldName,pXMLString->c_str(),i);
for(unsigned long j=i+1;j<nDims;j++) {
fprintf(outfile,"*_%s_lattice%li",fieldName,j);
}
fprintf(outfile,"*_%s_%s_ncomponents;\n",fieldName,pXMLString->c_str());
}
fprintf(outfile,"\n");
}
fprintf(outfile,"%s// copy cross vectors into K and I vectors\n",indent);
fprintf(outfile,"\n");
for(list<XMLString>::const_iterator pXMLString = crossVectorNamesList()->begin(); pXMLString != crossVectorNamesList()->end(); pXMLString++) {
fprintf(outfile,"%sfor(unsigned long _j=0; _j<_%s_cross_size*_%s_%s_ncomponents; _j++) {\n",
indent,fieldName,fieldName,pXMLString->c_str());
fprintf(outfile,"%s _%s_%s_K[_j] = _%s_%s[_%s_%s_index_pointer_begin + _j];\n",
indent,fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str());
fprintf(outfile,"%s _%s_%s_I[_j] = _%s_%s[_%s_%s_index_pointer_begin + _j];\n",
indent,fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str());
fprintf(outfile,"%s }\n",indent);
fprintf(outfile,"\n");
}
fprintf(outfile,"%s// store main vectors into old\n",indent);
fprintf(outfile,"\n");
for(list<XMLString>::const_iterator pXMLString = myMainVectorNamesList.begin(); pXMLString != myMainVectorNamesList.end(); pXMLString++) {
fprintf(outfile,"%sfor(unsigned long _j=0; _j<_%s_cross_size*_%s_%s_ncomponents; _j++)\n",
indent,fieldName,fieldName,pXMLString->c_str());
fprintf(outfile,"%s _%s_%s_old[_j] = _active_%s_%s[_%s_%s_index_pointer_begin + _j];\n",
indent,fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str());
fprintf(outfile,"\n");
}
// ****************************************************************
fprintf(outfile,"%s// ********** step 1 ***************\n",indent);
fprintf(outfile,"\n");
fprintf(outfile,"%s {\n",indent);
fprintf(outfile,"\n");
for(list<XMLString>::const_iterator pXMLString = myTotalVectorsList.begin(); pXMLString != myTotalVectorsList.end(); pXMLString++) {
fprintf(outfile,"%s unsigned long _%s_%s_index_pointer=_%s_%s_index_pointer_begin;\n",
indent,fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str());
}
fprintf(outfile,"\n");
fprintf(outfile,"%s // calculate k1\n",indent);
fprintf(outfile,"\n");
for(list<XMLString>::const_iterator pXMLString = crossVectorNamesList()->begin(); pXMLString != crossVectorNamesList()->end(); pXMLString++) {
fprintf(outfile,"%s unsigned long _%s_%s_index_pointer_local=0;\n",indent,fieldName,pXMLString->c_str());
}
fprintf(outfile,"\n");
// open inner loops
for(unsigned long i=crossDimNumber()+1; i<nDims; i++) {
for(unsigned long j=0; j<i; j++) {
fprintf(outfile," ");
}
fprintf(outfile," double %s = _%s_xmin%li;\n",simulation()->field()->geometry()->dimension(i)->name.c_str(),fieldName,i);
fprintf(outfile,"\n");
for(unsigned long j=0; j<i; j++) {
fprintf(outfile," ");
}
fprintf(outfile," for(unsigned long _i%li=0; _i%li<_%s_lattice%li; _i%li++) {\n",i,i,fieldName,i,i);
fprintf(outfile,"\n");
}
fprintf(outfile,"// *********** cross_propagation code ***********\n");
fprintf(outfile,"%s\n",crossPropagationCode()->c_str());
fprintf(outfile,"// **********************************************\n");
fprintf(outfile,"\n");
for(list<const xmdsVector*>::const_iterator pxmdsVector = crossVectorList.begin(); pxmdsVector != crossVectorList.end(); pxmdsVector++) {
for(unsigned long i=0;i<(*pxmdsVector)->nComponents();i++) {
fprintf(outfile,"%s _%s_%s_d[_%s_%s_index_pointer_local + %li] = d%s_d%s*_%s_dx%li;\n",
indent2,fieldName,(*pxmdsVector)->name()->c_str(),fieldName,(*pxmdsVector)->name()->c_str(),i,
(*pxmdsVector)->componentName(i)->c_str(),
simulation()->field()->geometry()->dimension(crossDimNumber())->name.c_str(),fieldName,crossDimNumber());
}
fprintf(outfile,"\n");
}
//close inner loops
if(crossDimNumber()+1<nDims) {
for(list<XMLString>::const_iterator pXMLString = myTotalVectorsList.begin(); pXMLString != myTotalVectorsList.end(); pXMLString++) {
fprintf(outfile,"%s _%s_%s_index_pointer += _%s_%s_ncomponents;\n",
indent2,fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str());
}
fprintf(outfile,"\n");
for(list<XMLString>::const_iterator pXMLString = crossVectorNamesList()->begin(); pXMLString != crossVectorNamesList()->end(); pXMLString++) {
fprintf(outfile,"%s _%s_%s_index_pointer_local += _%s_%s_ncomponents;\n",
indent2,fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str());
}
fprintf(outfile,"\n");
}
for(unsigned long i=nDims; i>crossDimNumber()+1; i--) {
for(unsigned long j=0; j<i; j++) {
fprintf(outfile," ");
}
fprintf(outfile," %s += _%s_dx%li;\n",simulation()->field()->geometry()->dimension(i-1)->name.c_str(),fieldName,i-1);
fprintf(outfile,"\n");
for(unsigned long j=0; j<i; j++) {
fprintf(outfile," ");
}
fprintf(outfile," }\n");
fprintf(outfile,"\n");
}
fprintf(outfile,"%s }\n",indent);
fprintf(outfile,"\n");
fprintf(outfile,"%s// add to K cross vector\n",indent);
fprintf(outfile,"\n");
for(list<XMLString>::const_iterator pXMLString = crossVectorNamesList()->begin(); pXMLString != crossVectorNamesList()->end(); pXMLString++) {
fprintf(outfile,"%sfor(unsigned long _j=0; _j<_%s_cross_size*_%s_%s_ncomponents; _j++)\n",
indent,fieldName,fieldName,pXMLString->c_str());
fprintf(outfile,"%s _%s_%s_K[_j] += _%s_%s_d[_j]/6;\n",indent,fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str());
fprintf(outfile,"\n");
}
fprintf(outfile,"%s// create next cross vectors\n",indent);
fprintf(outfile,"\n");
for(list<XMLString>::const_iterator pXMLString = crossVectorNamesList()->begin(); pXMLString != crossVectorNamesList()->end(); pXMLString++) {
fprintf(outfile,"%sfor(unsigned long _j=0; _j<_%s_cross_size*_%s_%s_ncomponents; _j++)\n",
indent,fieldName,fieldName,pXMLString->c_str());
fprintf(outfile,"%s _%s_%s[_%s_%s_index_pointer_begin + _j] = _%s_%s_I[_j] + _%s_%s_d[_j]/2;\n",
indent,fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str());
fprintf(outfile,"\n");
}
fprintf(outfile,"%s// create midpoint main vectors for steps 2 and 3\n",indent);
fprintf(outfile,"\n");
for(list<XMLString>::const_iterator pXMLString = myMainVectorNamesList.begin(); pXMLString != myMainVectorNamesList.end(); pXMLString++) {
fprintf(outfile,"%sfor(unsigned long _j=_%s_%s_index_pointer_begin; _j < _%s_%s_index_pointer_begin + _%s_cross_size*_%s_%s_ncomponents; _j++)\n",
indent,fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str(),fieldName,fieldName,pXMLString->c_str());
fprintf(outfile,"%s _active_%s_%s[_j] = (_active_%s_%s[_j]+_active_%s_%s[_%s_cross_size*_%s_%s_ncomponents + _j])/2;\n",
indent,fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str(),
fieldName,fieldName,pXMLString->c_str());
fprintf(outfile,"\n");
}
fprintf(outfile,"%s// move cross dim to lattice midpoint for steps 2 and 3\n",indent);
fprintf(outfile,"\n");
fprintf(outfile,"%s%s += _%s_dx%li/2;\n",
indent,simulation()->field()->geometry()->dimension(crossDimNumber())->name.c_str(),fieldName,crossDimNumber());
fprintf(outfile,"\n");
// ****************************************************************
fprintf(outfile,"%s// ********** step 2 ***************\n",indent);
fprintf(outfile,"\n");
fprintf(outfile,"%s {\n",indent);
fprintf(outfile,"\n");
for(list<XMLString>::const_iterator pXMLString = myTotalVectorsList.begin(); pXMLString != myTotalVectorsList.end(); pXMLString++) {
fprintf(outfile,"%s unsigned long _%s_%s_index_pointer=_%s_%s_index_pointer_begin;\n",
indent,fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str());
}
fprintf(outfile,"\n");
fprintf(outfile,"%s // calculate k1\n",indent);
fprintf(outfile,"\n");
for(list<XMLString>::const_iterator pXMLString = crossVectorNamesList()->begin(); pXMLString != crossVectorNamesList()->end(); pXMLString++) {
fprintf(outfile,"%s unsigned long _%s_%s_index_pointer_local=0;\n",indent,fieldName,pXMLString->c_str());
}
fprintf(outfile,"\n");
// open inner loops
for(unsigned long i=crossDimNumber()+1; i<nDims; i++) {
for(unsigned long j=0; j<i; j++) {
fprintf(outfile," ");
}
fprintf(outfile," double %s = _%s_xmin%li;\n",simulation()->field()->geometry()->dimension(i)->name.c_str(),fieldName,i);
fprintf(outfile,"\n");
for(unsigned long j=0; j<i; j++) {
fprintf(outfile," ");
}
fprintf(outfile," for(unsigned long _i%li=0; _i%li<_%s_lattice%li; _i%li++) {\n",i,i,fieldName,i,i);
fprintf(outfile,"\n");
}
fprintf(outfile,"// *********** cross_propagation code ***********\n");
fprintf(outfile,"%s\n",crossPropagationCode()->c_str());
fprintf(outfile,"// **********************************************\n");
fprintf(outfile,"\n");
for(list<const xmdsVector*>::const_iterator pxmdsVector = crossVectorList.begin(); pxmdsVector != crossVectorList.end(); pxmdsVector++) {
for(unsigned long i=0;i<(*pxmdsVector)->nComponents();i++) {
fprintf(outfile,"%s _%s_%s_d[_%s_%s_index_pointer_local + %li] = d%s_d%s*_%s_dx%li;\n",
indent2,fieldName,(*pxmdsVector)->name()->c_str(),fieldName,(*pxmdsVector)->name()->c_str(),i,
(*pxmdsVector)->componentName(i)->c_str(),
simulation()->field()->geometry()->dimension(crossDimNumber())->name.c_str(),fieldName,crossDimNumber());
}
fprintf(outfile,"\n");
}
//close inner loops
if(crossDimNumber()+1<nDims) {
for(list<XMLString>::const_iterator pXMLString = myTotalVectorsList.begin(); pXMLString != myTotalVectorsList.end(); pXMLString++) {
fprintf(outfile,"%s _%s_%s_index_pointer += _%s_%s_ncomponents;\n",
indent2,fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str());
}
fprintf(outfile,"\n");
for(list<XMLString>::const_iterator pXMLString = crossVectorNamesList()->begin(); pXMLString != crossVectorNamesList()->end(); pXMLString++) {
fprintf(outfile,"%s _%s_%s_index_pointer_local += _%s_%s_ncomponents;\n",
indent2,fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str());
}
fprintf(outfile,"\n");
}
for(unsigned long i=nDims; i>crossDimNumber()+1; i--) {
for(unsigned long j=0; j<i; j++) {
fprintf(outfile," ");
}
fprintf(outfile," %s += _%s_dx%li;\n",simulation()->field()->geometry()->dimension(i-1)->name.c_str(),fieldName,i-1);
fprintf(outfile,"\n");
for(unsigned long j=0; j<i; j++) {
fprintf(outfile," ");
}
fprintf(outfile," }\n");
fprintf(outfile,"\n");
}
fprintf(outfile,"%s }\n",indent);
fprintf(outfile,"\n");
fprintf(outfile,"%s// add to K cross vector\n",indent);
fprintf(outfile,"\n");
for(list<XMLString>::const_iterator pXMLString = crossVectorNamesList()->begin(); pXMLString != crossVectorNamesList()->end(); pXMLString++) {
fprintf(outfile,"%sfor(unsigned long _j=0; _j<_%s_cross_size*_%s_%s_ncomponents; _j++)\n",
indent,fieldName,fieldName,pXMLString->c_str());
fprintf(outfile,"%s _%s_%s_K[_j] += _%s_%s_d[_j]/3;\n",indent,fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str());
fprintf(outfile,"\n");
}
fprintf(outfile,"%s// create next cross vectors\n",indent);
fprintf(outfile,"\n");
for(list<XMLString>::const_iterator pXMLString = crossVectorNamesList()->begin(); pXMLString != crossVectorNamesList()->end(); pXMLString++) {
fprintf(outfile,"%sfor(unsigned long _j=0; _j<_%s_cross_size*_%s_%s_ncomponents; _j++)\n",
indent,fieldName,fieldName,pXMLString->c_str());
fprintf(outfile,"%s _%s_%s[_%s_%s_index_pointer_begin + _j] = _%s_%s_I[_j] + _%s_%s_d[_j]/2;\n",
indent,fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str());
fprintf(outfile,"\n");
}
// ****************************************************************
fprintf(outfile,"%s// ********** step 3 ***************\n",indent);
fprintf(outfile,"\n");
fprintf(outfile,"%s {\n",indent);
fprintf(outfile,"\n");
for(list<XMLString>::const_iterator pXMLString = myTotalVectorsList.begin(); pXMLString != myTotalVectorsList.end(); pXMLString++) {
fprintf(outfile,"%s unsigned long _%s_%s_index_pointer=_%s_%s_index_pointer_begin;\n",
indent,fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str());
}
fprintf(outfile,"\n");
fprintf(outfile,"%s // calculate k1\n",indent);
fprintf(outfile,"\n");
for(list<XMLString>::const_iterator pXMLString = crossVectorNamesList()->begin(); pXMLString != crossVectorNamesList()->end(); pXMLString++) {
fprintf(outfile,"%s unsigned long _%s_%s_index_pointer_local=0;\n",indent,fieldName,pXMLString->c_str());
}
fprintf(outfile,"\n");
// open inner loops
for(unsigned long i=crossDimNumber()+1; i<nDims; i++) {
for(unsigned long j=0; j<i; j++) {
fprintf(outfile," ");
}
fprintf(outfile," double %s = _%s_xmin%li;\n",simulation()->field()->geometry()->dimension(i)->name.c_str(),fieldName,i);
fprintf(outfile,"\n");
for(unsigned long j=0; j<i; j++) {
fprintf(outfile," ");
}
fprintf(outfile," for(unsigned long _i%li=0; _i%li<_%s_lattice%li; _i%li++) {\n",i,i,fieldName,i,i);
fprintf(outfile,"\n");
}
fprintf(outfile,"// *********** cross_propagation code ***********\n");
fprintf(outfile,"%s\n",crossPropagationCode()->c_str());
fprintf(outfile,"// **********************************************\n");
fprintf(outfile,"\n");
for(list<const xmdsVector*>::const_iterator pxmdsVector = crossVectorList.begin(); pxmdsVector != crossVectorList.end(); pxmdsVector++) {
for(unsigned long i=0;i<(*pxmdsVector)->nComponents();i++) {
fprintf(outfile,"%s _%s_%s_d[_%s_%s_index_pointer_local + %li] = d%s_d%s*_%s_dx%li;\n",
indent2,fieldName,(*pxmdsVector)->name()->c_str(),fieldName,(*pxmdsVector)->name()->c_str(),i,
(*pxmdsVector)->componentName(i)->c_str(),
simulation()->field()->geometry()->dimension(crossDimNumber())->name.c_str(),fieldName,crossDimNumber());
}
fprintf(outfile,"\n");
}
//close inner loops
if(crossDimNumber()+1<nDims) {
for(list<XMLString>::const_iterator pXMLString = myTotalVectorsList.begin(); pXMLString != myTotalVectorsList.end(); pXMLString++) {
fprintf(outfile,"%s _%s_%s_index_pointer += _%s_%s_ncomponents;\n",
indent2,fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str());
}
fprintf(outfile,"\n");
for(list<XMLString>::const_iterator pXMLString = crossVectorNamesList()->begin(); pXMLString != crossVectorNamesList()->end(); pXMLString++) {
fprintf(outfile,"%s _%s_%s_index_pointer_local += _%s_%s_ncomponents;\n",
indent2,fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str());
}
fprintf(outfile,"\n");
}
for(unsigned long i=nDims; i>crossDimNumber()+1; i--) {
for(unsigned long j=0; j<i; j++) {
fprintf(outfile," ");
}
fprintf(outfile," %s += _%s_dx%li;\n",simulation()->field()->geometry()->dimension(i-1)->name.c_str(),fieldName,i-1);
fprintf(outfile,"\n");
for(unsigned long j=0; j<i; j++) {
fprintf(outfile," ");
}
fprintf(outfile," }\n");
fprintf(outfile,"\n");
}
fprintf(outfile,"%s }\n",indent);
fprintf(outfile,"\n");
fprintf(outfile,"%s// add to K cross vector\n",indent);
fprintf(outfile,"\n");
for(list<XMLString>::const_iterator pXMLString = crossVectorNamesList()->begin(); pXMLString != crossVectorNamesList()->end(); pXMLString++) {
fprintf(outfile,"%sfor(unsigned long _j=0; _j<_%s_cross_size*_%s_%s_ncomponents; _j++)\n",
indent,fieldName,fieldName,pXMLString->c_str());
fprintf(outfile,"%s _%s_%s_K[_j] += _%s_%s_d[_j]/3;\n",indent,fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str());
fprintf(outfile,"\n");
}
fprintf(outfile,"%s// create next cross vectors\n",indent);
fprintf(outfile,"\n");
for(list<XMLString>::const_iterator pXMLString = crossVectorNamesList()->begin(); pXMLString != crossVectorNamesList()->end(); pXMLString++) {
fprintf(outfile,"%sfor(unsigned long _j=0; _j<_%s_cross_size*_%s_%s_ncomponents; _j++)\n",
indent,fieldName,fieldName,pXMLString->c_str());
fprintf(outfile,"%s _%s_%s[_%s_%s_index_pointer_begin + _j] = _%s_%s_I[_j] + _%s_%s_d[_j];\n",
indent,fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str());
fprintf(outfile,"\n");
}
fprintf(outfile,"%s// move main vectors to next lattice point for step 4\n",indent);
fprintf(outfile,"\n");
for(list<XMLString>::const_iterator pXMLString = myMainVectorNamesList.begin(); pXMLString != myMainVectorNamesList.end(); pXMLString++) {
fprintf(outfile,"%sfor(unsigned long _j=_%s_%s_index_pointer_begin; _j < _%s_%s_index_pointer_begin + _%s_cross_size*_%s_%s_ncomponents; _j++)\n",
indent,fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str(),fieldName,fieldName,pXMLString->c_str());
fprintf(outfile,"%s _active_%s_%s[_j] = _active_%s_%s[_%s_cross_size*_%s_%s_ncomponents + _j];\n",
indent,fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str(),
fieldName,fieldName,pXMLString->c_str());
fprintf(outfile,"\n");
}
fprintf(outfile,"%s// move cross dim to next lattice point for step 4\n",indent);
fprintf(outfile,"\n");
fprintf(outfile,"%s%s += _%s_dx%li/2;\n",
indent,simulation()->field()->geometry()->dimension(crossDimNumber())->name.c_str(),fieldName,crossDimNumber());
fprintf(outfile,"\n");
// ****************************************************************
fprintf(outfile,"%s// ********** step 4 ***************\n",indent);
fprintf(outfile,"\n");
fprintf(outfile,"%s {\n",indent);
fprintf(outfile,"\n");
for(list<XMLString>::const_iterator pXMLString = myTotalVectorsList.begin(); pXMLString != myTotalVectorsList.end(); pXMLString++) {
fprintf(outfile,"%s unsigned long _%s_%s_index_pointer=_%s_%s_index_pointer_begin;\n",
indent,fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str());
}
fprintf(outfile,"\n");
fprintf(outfile,"%s // calculate k1\n",indent);
fprintf(outfile,"\n");
for(list<XMLString>::const_iterator pXMLString = crossVectorNamesList()->begin(); pXMLString != crossVectorNamesList()->end(); pXMLString++) {
fprintf(outfile,"%s unsigned long _%s_%s_index_pointer_local=0;\n",indent,fieldName,pXMLString->c_str());
}
fprintf(outfile,"\n");
// open inner loops
for(unsigned long i=crossDimNumber()+1; i<nDims; i++) {
for(unsigned long j=0; j<i; j++) {
fprintf(outfile," ");
}
fprintf(outfile," double %s = _%s_xmin%li;\n",simulation()->field()->geometry()->dimension(i)->name.c_str(),fieldName,i);
fprintf(outfile,"\n");
for(unsigned long j=0; j<i; j++) {
fprintf(outfile," ");
}
fprintf(outfile," for(unsigned long _i%li=0; _i%li<_%s_lattice%li; _i%li++) {\n",i,i,fieldName,i,i);
fprintf(outfile,"\n");
}
fprintf(outfile,"// *********** cross_propagation code ***********\n");
fprintf(outfile,"%s\n",crossPropagationCode()->c_str());
fprintf(outfile,"// **********************************************\n");
fprintf(outfile,"\n");
for(list<const xmdsVector*>::const_iterator pxmdsVector = crossVectorList.begin(); pxmdsVector != crossVectorList.end(); pxmdsVector++) {
for(unsigned long i=0;i<(*pxmdsVector)->nComponents();i++) {
fprintf(outfile,"%s _%s_%s_d[_%s_%s_index_pointer_local + %li] = d%s_d%s*_%s_dx%li;\n",
indent2,fieldName,(*pxmdsVector)->name()->c_str(),fieldName,(*pxmdsVector)->name()->c_str(),i,
(*pxmdsVector)->componentName(i)->c_str(),
simulation()->field()->geometry()->dimension(crossDimNumber())->name.c_str(),fieldName,crossDimNumber());
}
fprintf(outfile,"\n");
}
//close inner loops
if(crossDimNumber()+1<nDims) {
for(list<XMLString>::const_iterator pXMLString = myTotalVectorsList.begin(); pXMLString != myTotalVectorsList.end(); pXMLString++) {
fprintf(outfile,"%s _%s_%s_index_pointer += _%s_%s_ncomponents;\n",
indent2,fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str());
}
fprintf(outfile,"\n");
for(list<XMLString>::const_iterator pXMLString = crossVectorNamesList()->begin(); pXMLString != crossVectorNamesList()->end(); pXMLString++) {
fprintf(outfile,"%s _%s_%s_index_pointer_local += _%s_%s_ncomponents;\n",
indent2,fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str());
}
fprintf(outfile,"\n");
}
for(unsigned long i=nDims; i>crossDimNumber()+1; i--) {
for(unsigned long j=0; j<i; j++) {
fprintf(outfile," ");
}
fprintf(outfile," %s += _%s_dx%li;\n",simulation()->field()->geometry()->dimension(i-1)->name.c_str(),fieldName,i-1);
fprintf(outfile,"\n");
for(unsigned long j=0; j<i; j++) {
fprintf(outfile," ");
}
fprintf(outfile," }\n");
fprintf(outfile,"\n");
}
fprintf(outfile,"%s }\n",indent);
fprintf(outfile,"\n");
fprintf(outfile,"%s// add to K cross vector\n",indent);
fprintf(outfile,"\n");
for(list<XMLString>::const_iterator pXMLString = crossVectorNamesList()->begin(); pXMLString != crossVectorNamesList()->end(); pXMLString++) {
fprintf(outfile,"%sfor(unsigned long _j=0; _j<_%s_cross_size*_%s_%s_ncomponents; _j++)\n",
indent,fieldName,fieldName,pXMLString->c_str());
fprintf(outfile,"%s _%s_%s_K[_j] += _%s_%s_d[_j]/6;\n",indent,fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str());
fprintf(outfile,"\n");
}
fprintf(outfile,"%s// copy I cross vector back into old main cross vector\n",indent);
fprintf(outfile,"\n");
for(list<XMLString>::const_iterator pXMLString = crossVectorNamesList()->begin(); pXMLString != crossVectorNamesList()->end(); pXMLString++) {
fprintf(outfile,"%sfor(unsigned long _j=0; _j<_%s_cross_size*_%s_%s_ncomponents; _j++)\n",
indent,fieldName,fieldName,pXMLString->c_str());
fprintf(outfile,"%s _%s_%s[_%s_%s_index_pointer_begin + _j] = _%s_%s_I[_j];\n",
indent,fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str());
fprintf(outfile,"\n");
}
fprintf(outfile,"%s// copy K cross vector into next main cross vector\n",indent);
fprintf(outfile,"\n");
for(list<XMLString>::const_iterator pXMLString = crossVectorNamesList()->begin(); pXMLString != crossVectorNamesList()->end(); pXMLString++) {
fprintf(outfile,"%sfor(unsigned long _j=0; _j<_%s_cross_size*_%s_%s_ncomponents; _j++)\n",
indent,fieldName,fieldName,pXMLString->c_str());
fprintf(outfile,"%s _%s_%s[_%s_%s_index_pointer_begin + _%s_cross_size*_%s_%s_ncomponents + _j] = _%s_%s_K[_j];\n",
indent,fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str(),fieldName,fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str());
fprintf(outfile,"\n");
}
fprintf(outfile,"%s// copy old main vectors back into last lattice point\n",indent);
fprintf(outfile,"\n");
for(list<XMLString>::const_iterator pXMLString = myMainVectorNamesList.begin(); pXMLString != myMainVectorNamesList.end(); pXMLString++) {
fprintf(outfile,"%sfor(unsigned long _j=0; _j<_%s_cross_size*_%s_%s_ncomponents; _j++)\n",
indent,fieldName,fieldName,pXMLString->c_str());
fprintf(outfile,"%s _active_%s_%s[_%s_%s_index_pointer_begin + _j] = _%s_%s_old[_j];\n",
indent,fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str());
fprintf(outfile,"\n");
}
// close outer loops
fprintf(outfile,"%s}\n",indent);
for(unsigned long i=crossDimNumber(); i>0; i--) {
fprintf(outfile,"\n");
for(unsigned long j=0; j<i; j++) {
fprintf(outfile," ");
}
fprintf(outfile,"%s += _%s_dx%li;\n",simulation()->field()->geometry()->dimension(i-1)->name.c_str(),fieldName,i-1);
fprintf(outfile,"\n");
for(unsigned long j=0; j<i; j++) {
fprintf(outfile," ");
}
fprintf(outfile,"}\n");
}
for(list<XMLString>::const_iterator pXMLString = myMainVectorNamesList.begin(); pXMLString != myMainVectorNamesList.end(); pXMLString++) {
fprintf(outfile," delete[] _%s_%s_old;\n", fieldName, pXMLString->c_str());
fprintf(outfile,"\n");
}
for (list<XMLString>::const_iterator pXMLString = crossVectorNamesList()->begin(); pXMLString != crossVectorNamesList()->end(); pXMLString++) {
fprintf(outfile, "delete[] _%s_%s_K;\n", fieldName, pXMLString->c_str());
fprintf(outfile, "delete[] _%s_%s_I;\n", fieldName, pXMLString->c_str());
fprintf(outfile, "delete[] _%s_%s_d;\n", fieldName, pXMLString->c_str());
fprintf(outfile,"\n");
}
fprintf(outfile,"\n");
fprintf(outfile,"}\n");
fprintf(outfile,"\n");
};
syntax highlighted by Code2HTML, v. 0.9.1