xmdsintegratesiip.cc

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

/*
  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.24 2004/10/06 07:52:38 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");
  }