xmdsintegraterk4.cc

XMDS is a code generator that integrates equations. You write them down in human readable form in a

    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");
};