/*
 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: xmdsintegraterk4ex.cc,v 1.13 2004/09/15 06:01:58 joehope Exp $
*/

/*! @file xmdsintegraterk4ex.cc
  @brief Integrate element parsing classes and methods; fourth order Runge-Kutta in the explicit picture

  More detailed explanation...
*/

#include<xmlbasics.h>
#include<dom3.h>
#include<xmdsutils.h>
#include<xmdsclasses.h>

// ******************************************************************************
// ******************************************************************************
//                              xmdsIntegrateRK4EX public
// ******************************************************************************
// ******************************************************************************

extern bool debugFlag;

long nxmdsIntegrateRK4EXs=0;  //!< The number of integrate RK4EX objects

// ******************************************************************************
xmdsIntegrateRK4EX::xmdsIntegrateRK4EX(
				       const xmdsSimulation *const yourSimulation,
				       const bool& yourVerboseMode) :
  xmdsIntegrate(yourSimulation,yourVerboseMode),
  xmdsIntegrateEX(yourSimulation,yourVerboseMode),
  xmdsIntegrateRK4(yourSimulation,yourVerboseMode) {
  if(debugFlag) {
    nxmdsIntegrateRK4EXs++;
    printf("xmdsIntegrateRK4EX::xmdsIntegrateRK4EX\n");
    printf("nxmdsIntegrateRK4EXs=%li\n",nxmdsIntegrateRK4EXs);
  }
};

// ******************************************************************************
xmdsIntegrateRK4EX::~xmdsIntegrateRK4EX() {
  if(debugFlag) {
    nxmdsIntegrateRK4EXs--;
    printf("xmdsIntegrateRK4EX::~xmdsIntegrateRK4EX\n");
    printf("nxmdsIntegrateRK4EXs=%li\n",nxmdsIntegrateRK4EXs);
  }
};

// ******************************************************************************
void xmdsIntegrateRK4EX::processElement(
					const Element *const yourElement) {
  if(debugFlag) {
    printf("xmdsIntegrateRK4EX::processElement\n");
  }

  if(verbose()) {
    printf("Processing integrate RK4EX element ...\n");
  }

  xmdsIntegrate::processElement(yourElement);
  xmdsIntegrateRK4::processElement(yourElement);
  xmdsIntegrateEX::processElement(yourElement);
};

// ******************************************************************************
// ******************************************************************************
//                              xmdsIntegrateRK4EX private
// ******************************************************************************
// ******************************************************************************

// ******************************************************************************
void xmdsIntegrateRK4EX::writePrototypes(
					 FILE *const outfile) const {
  if(debugFlag) {
    printf("xmdsIntegrateRK4EX::writePrototypes\n");
  }

  fprintf(outfile,"// ********************************************************\n");
  fprintf(outfile,"// segment %li (RK4EX) prototypes\n",segmentNumber);
  fprintf(outfile,"\n");

  xmdsIntegrate::writePrototypes(outfile);
  xmdsIntegrateEX::writePrototypes(outfile);
  xmdsIntegrateRK4::writePrototypes(outfile);
};

// ******************************************************************************
void xmdsIntegrateRK4EX::writeRoutines(
				       FILE *const outfile) const {
  if(debugFlag) {
    printf("xmdsIntegrateRK4EX::writeRoutines\n");
  }

  fprintf(outfile,"// ********************************************************\n");
  fprintf(outfile,"// segment %li (RK4EX) routines\n",segmentNumber);
  fprintf(outfile,"\n");

  xmdsIntegrate::writeRoutines(outfile);
  xmdsIntegrateEX::writeRoutines(outfile);
  xmdsIntegrateRK4::writeRoutines(outfile);
};

// ******************************************************************************
void xmdsIntegrateRK4EX::writeSingleStepCode(
					     FILE *const outfile,
					     const stepCaseEnum& stepCase) const {
  if(debugFlag) {
    printf("xmdsIntegrateRK4EX::writeSingleStepCode\n");
  }

  const char *const fieldName = simulation()->field()->name()->c_str();
  const char *const propDim = simulation()->parameters()->propDimName.c_str();

  list<XMLString> tempVectorNamesList;
  tempVectorNamesList.push_back("main");

  const char* indent = "	";
  if(simulation()->parameters()->errorCheck) {
    indent = "		";
  }

  const char* noiseVector="";
  if((simulation()->parameters()->stochastic)&&(!noNoises())) {
    noiseVector=",_noise_vector";
  }

  simulation()->field()->vectors2space(outfile,0,tempVectorNamesList,indent);

  fprintf(outfile,"%s// a_k=a\n",indent);
  if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
    fprintf(outfile,"%sfor(long _i1=0; _i1<total_local_size*_%s_main_ncomponents; _i1++)\n",indent,fieldName);
  }
  else {
    fprintf(outfile,"%sfor(unsigned long _i1=0; _i1<_%s_size*_%s_main_ncomponents; _i1++)\n",indent,fieldName,fieldName);
  }
  fprintf(outfile,"%s	akfield_%s_main[_i1] = _%s_main[_i1];\n",indent,fieldName,fieldName);
  fprintf(outfile,"\n");

  fprintf(outfile,"%s// a_i=a\n",indent);
  if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
    fprintf(outfile,"%sfor(long _i1=0; _i1<total_local_size*_%s_main_ncomponents; _i1++)\n",indent,fieldName);
  }
  else {
    fprintf(outfile,"%sfor(unsigned long _i1=0; _i1<_%s_size*_%s_main_ncomponents; _i1++)\n",indent,fieldName,fieldName);
  }
  fprintf(outfile,"%s	aifield_%s_main[_i1] = _%s_main[_i1];\n",indent,fieldName,fieldName);
  fprintf(outfile,"\n");

  fprintf(outfile,"%s_active_%s_main = akfield_%s_main;\n",indent,fieldName,fieldName);
  fprintf(outfile,"\n");

  if(usesKOperators()) {
    fprintf(outfile,"%s// calculate the co derivative terms\n",indent);
    fprintf(outfile,"%s_segment%li_calculate_co_terms();\n",indent,segmentNumber);
    fprintf(outfile,"\n");
  }

  fprintf(outfile,"%s// a_k=G[a_k,t]\n",indent);
  fprintf(outfile,"%s_segment%li_calculate_delta_a(_step,cycle%s);\n",indent,segmentNumber,noiseVector);
  fprintf(outfile,"\n");

  fprintf(outfile,"%s// a=a+a_k/6\n",indent);
  if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
    fprintf(outfile,"%sfor(long _i1=0; _i1<total_local_size*_%s_main_ncomponents; _i1++)\n",indent,fieldName);
  }
  else {
    fprintf(outfile,"%sfor(unsigned long _i1=0; _i1<_%s_size*_%s_main_ncomponents; _i1++)\n",indent,fieldName,fieldName);
  }
  fprintf(outfile,"%s	_%s_main[_i1] += akfield_%s_main[_i1]/6;\n",indent,fieldName,fieldName);
  fprintf(outfile,"\n");

  fprintf(outfile,"%s// a_k = a_i + a_k/2\n",indent);
  if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
    fprintf(outfile,"%sfor(long _i1=0; _i1<total_local_size*_%s_main_ncomponents; _i1++)\n",indent,fieldName);
  }
  else {
    fprintf(outfile,"%sfor(unsigned long _i1=0; _i1<_%s_size*_%s_main_ncomponents; _i1++)\n",indent,fieldName,fieldName);
  }
  fprintf(outfile,"%s	akfield_%s_main[_i1] = aifield_%s_main[_i1] + akfield_%s_main[_i1]/2;\n",indent,fieldName,fieldName,fieldName);
  fprintf(outfile,"\n");

  fprintf(outfile,"%s%s += _step/2;\n",indent,propDim);
  fprintf(outfile,"\n");

  if(usesKOperators()) {
    fprintf(outfile,"%s// re-calculate the co derivative terms\n",indent);
    fprintf(outfile,"%s_segment%li_calculate_co_terms();\n",indent,segmentNumber);
    fprintf(outfile,"\n");
  }

  fprintf(outfile,"%s// a_k=G[a_k,t+h/2]\n",indent);
  fprintf(outfile,"%s_segment%li_calculate_delta_a(_step,cycle%s);\n",indent,segmentNumber,noiseVector);
  fprintf(outfile,"\n");

  fprintf(outfile,"%s// a=a+a_k/3\n",indent);
  if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
    fprintf(outfile,"%sfor(long _i1=0; _i1<total_local_size*_%s_main_ncomponents; _i1++)\n",indent,fieldName);
  }
  else {
    fprintf(outfile,"%sfor(unsigned long _i1=0; _i1<_%s_size*_%s_main_ncomponents; _i1++)\n",indent,fieldName,fieldName);
  }
  fprintf(outfile,"%s	_%s_main[_i1] += akfield_%s_main[_i1]/3;\n",indent,fieldName,fieldName);
  fprintf(outfile,"\n");

  fprintf(outfile,"%s// a_k = a_i + a_k/2\n",indent);
  if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
    fprintf(outfile,"%sfor(long _i1=0; _i1<total_local_size*_%s_main_ncomponents; _i1++)\n",indent,fieldName);
  }
  else {
    fprintf(outfile,"%sfor(unsigned long _i1=0; _i1<_%s_size*_%s_main_ncomponents; _i1++)\n",indent,fieldName,fieldName);
  }
  fprintf(outfile,"%s	akfield_%s_main[_i1] = aifield_%s_main[_i1] + akfield_%s_main[_i1]/2;\n",indent,fieldName,fieldName,fieldName);
  fprintf(outfile,"\n");

  if(usesKOperators()) {
    fprintf(outfile,"%s// re-calculate the co derivative terms\n",indent);
    fprintf(outfile,"%s_segment%li_calculate_co_terms();\n",indent,segmentNumber);
    fprintf(outfile,"\n");
  }

  fprintf(outfile,"%s// a_k=G[a_k,t+h/2]\n",indent);
  fprintf(outfile,"%s_segment%li_calculate_delta_a(_step,cycle%s);\n",indent,segmentNumber,noiseVector);
  fprintf(outfile,"\n");

  fprintf(outfile,"%s// a=a+a_k/3\n",indent);
  if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
    fprintf(outfile,"%sfor(long _i1=0; _i1<total_local_size*_%s_main_ncomponents; _i1++)\n",indent,fieldName);
  }
  else {
    fprintf(outfile,"%sfor(unsigned long _i1=0; _i1<_%s_size*_%s_main_ncomponents; _i1++)\n",indent,fieldName,fieldName);
  }
  fprintf(outfile,"%s	_%s_main[_i1] += akfield_%s_main[_i1]/3;\n",indent,fieldName,fieldName);
  fprintf(outfile,"\n");

  fprintf(outfile,"%s// a_k = a_i + a_k\n",indent);
  if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
    fprintf(outfile,"%sfor(long _i1=0; _i1<total_local_size*_%s_main_ncomponents; _i1++)\n",indent,fieldName);
  }
  else {
    fprintf(outfile,"%sfor(unsigned long _i1=0; _i1<_%s_size*_%s_main_ncomponents; _i1++)\n",indent,fieldName,fieldName);
  }
  fprintf(outfile,"%s	akfield_%s_main[_i1] = aifield_%s_main[_i1] + akfield_%s_main[_i1];\n",indent,fieldName,fieldName,fieldName);
  fprintf(outfile,"\n");

  fprintf(outfile,"%s%s += _step/2;\n",indent,propDim);
  fprintf(outfile,"\n");

  if(usesKOperators()) {
    fprintf(outfile,"%s// re-calculate the co derivative terms\n",indent);
    fprintf(outfile,"%s_segment%li_calculate_co_terms();\n",indent,segmentNumber);
    fprintf(outfile,"\n");
  }

  fprintf(outfile,"%s// a_k=G[a_k,t+h]\n",indent);
  fprintf(outfile,"%s_segment%li_calculate_delta_a(_step,cycle%s);\n",indent,segmentNumber,noiseVector);
  fprintf(outfile,"\n");

  fprintf(outfile,"%s_active_%s_main = _%s_main;\n",indent,fieldName,fieldName);
  fprintf(outfile,"\n");

  fprintf(outfile,"%s// a=a+a_k/6\n",indent);
  if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
    fprintf(outfile,"%sfor(long _i1=0; _i1<total_local_size*_%s_main_ncomponents; _i1++)\n",indent,fieldName);
  }
  else {
    fprintf(outfile,"%sfor(unsigned long _i1=0; _i1<_%s_size*_%s_main_ncomponents; _i1++)\n",indent,fieldName,fieldName);
  }
  fprintf(outfile,"%s	_%s_main[_i1] += akfield_%s_main[_i1]/6;\n",indent,fieldName,fieldName);
};