/*
 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: xmdsintegrateark45ex.cc,v 1.3 2005/05/02 02:35:58 sebwuester Exp $
*/

/*! @file xmdsintegrateark45ex.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>

// ******************************************************************************
// ******************************************************************************
//                              xmdsIntegrateARK45EX public
// ******************************************************************************
// ******************************************************************************

extern bool debugFlag;

long nxmdsIntegrateARK45EXs=0;  //!< The number of integrate ARK45EX objects

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

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

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

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

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

// ******************************************************************************
// ******************************************************************************
//                              xmdsIntegrateARK45EX private
// ******************************************************************************
// ******************************************************************************

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

  if(!simulation()->field()->getVector("main",mainVector)) {
    throw xmdsException("Internal error in xmdsIntegrateARKIP::writeResetRoutine: 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,"// segment %li (ARK45EX) prototypes\n",segmentNumber);
  fprintf(outfile,"\n");    
  xmdsIntegrate::writePrototypes(outfile);
  xmdsIntegrateEX::writePrototypes(outfile);
  
  fprintf(outfile,"\n");
  fprintf(outfile,"// integrate (ARK45EX) prototypes\n",segmentNumber);
  fprintf(outfile,"\n");
  fprintf(outfile,"void _segment%li_reset(%s* _reset_to,double _step);\n",segmentNumber,typeName);  
  xmdsIntegrateARK45::writePrototypes(outfile);
};

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

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

  xmdsIntegrate::writeRoutines(outfile);
  writeResetRoutine(outfile);
  xmdsIntegrateEX::writeRoutines(outfile);
  xmdsIntegrateARK45::writeRoutines(outfile);
};

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

  const char *const fieldName = simulation()->field()->name()->c_str();	
  
  const xmdsVector* mainVector;

  if(!simulation()->field()->getVector("main",mainVector)) {
    throw xmdsException("Internal error in xmdsIntegrateARK45::writeResetRoutine: 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_reset(%s* _reset_to,double _step){\n",segmentNumber,typeName);     
  if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
    fprintf(outfile,"for(long _i1=0; _i1<total_local_size*_%s_main_ncomponents; _i1++)\n",fieldName);
  }
  else {
    fprintf(outfile,"for(unsigned long _i1=0; _i1<_%s_size*_%s_main_ncomponents; _i1++)\n",fieldName,fieldName);
  }
  fprintf(outfile,"    _%s_main[_i1] = _reset_to[_i1];\n",fieldName);
  fprintf(outfile,"}\n");	
  fprintf(outfile,"\n");				     
					     
}

// ******************************************************************************
void xmdsIntegrateARK45EX::writeSingleStepCode(
					     FILE *const outfile,
					     const stepCaseEnum& stepCase) const {
  if(debugFlag) {
    printf("xmdsIntegrateARK45EX::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=y2=y1\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_check[_i1] = _%s_main[_i1];\n",indent,fieldName,fieldName,fieldName);
  fprintf(outfile,"\n");

  fprintf(outfile,"%s// a_i=y1\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// y1=y1+c_1*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	_%s_main[_i1] += c[1]*akfield_%s_main[_i1];\n",indent,fieldName,fieldName);
  
  fprintf(outfile,"%s// y2=y2+cs_1*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	_%s_check[_i1] += cs[1]*akfield_%s_main[_i1];\n",indent,fieldName,fieldName);

  fprintf(outfile,"%s// a_k = a_i + b_21*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] + b[2][1]*akfield_%s_main[_i1];\n",indent,fieldName,fieldName,fieldName);
  fprintf(outfile,"\n");

  fprintf(outfile,"%s%s += a[2]*_step;\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+aa_2*dt]\n",indent);
  fprintf(outfile,"%s_segment%li_calculate_delta_a(_step,cycle%s);\n",indent,segmentNumber,noiseVector);
  fprintf(outfile,"\n");
  
  fprintf(outfile,"%s// c_2==cs_2==0\n",indent);
  
    fprintf(outfile,"%s// a_j = d_1*a_i + d_2*y1 + d_3*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	ajfield_%s_main[_i1] = d[1]*aifield_%s_main[_i1] + d[2]*_%s_main[_i1]+d[3]*akfield_%s_main[_i1] ;\n",indent,fieldName,fieldName,fieldName,fieldName);
  fprintf(outfile,"\n");
  
  fprintf(outfile,"%s%s += a[3]*_step;\n",indent,propDim);
  fprintf(outfile,"\n");
  
  fprintf(outfile,"%s_active_%s_main = ajfield_%s_main;\n",indent,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_j=G[a_j,t+aa_3*dt]\n",indent);
  fprintf(outfile,"%s_segment%li_calculate_delta_a(_step,cycle%s);\n",indent,segmentNumber,noiseVector);
  fprintf(outfile,"\n");

  fprintf(outfile,"%s// a_l = e_1*a_i + e_2*y1 + e_3*a_k + e_4*a_j\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	alfield_%s_main[_i1] = e[1]*aifield_%s_main[_i1] + e[2]*_%s_main[_i1] + e[3]*akfield_%s_main[_i1] + e[4]*ajfield_%s_main[_i1] ;\n",indent,fieldName,fieldName,fieldName,fieldName,fieldName);
  fprintf(outfile,"\n");

  fprintf(outfile,"%s// y1=y1+c_3*a_j\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] += c[3]*ajfield_%s_main[_i1];\n",indent,fieldName,fieldName);
  
  fprintf(outfile,"%s// y2=y2+cs_3*a_j\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_check[_i1] += cs[3]*ajfield_%s_main[_i1];\n",indent,fieldName,fieldName);
  
  fprintf(outfile,"%s%s += a[4]*_step;\n",indent,propDim);
  fprintf(outfile,"\n");
  
  fprintf(outfile,"%s_active_%s_main = alfield_%s_main;\n",indent,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_l=G[a_l,t+aa_4*dt]\n",indent);
  fprintf(outfile,"%s_segment%li_calculate_delta_a(_step,cycle%s);\n",indent,segmentNumber,noiseVector);
  fprintf(outfile,"\n");

  fprintf(outfile,"%s// y1=y1+c_4*a_j\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] += c[4]*alfield_%s_main[_i1];\n",indent,fieldName,fieldName);
  
  fprintf(outfile,"%s// y2=y2+cs_4*a_j\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_check[_i1] += cs[4]*alfield_%s_main[_i1];\n",indent,fieldName,fieldName);

  fprintf(outfile,"%s// a_l = f_1*a_i + f_2*y1 + f_3*a_k + f_4*a_j +f_5*a_l;\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	alfield_%s_main[_i1] = f[1]*aifield_%s_main[_i1] + f[2]*_%s_main[_i1] + f[3]*akfield_%s_main[_i1] + f[4]*ajfield_%s_main[_i1] + f[5]*alfield_%s_main[_i1] ;\n",indent,fieldName,fieldName,fieldName,fieldName,fieldName,fieldName);
  fprintf(outfile,"\n");

  fprintf(outfile,"%s%s += a[5]*_step;\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_l=G[a_l,t+aa_5*dt]\n",indent);
  fprintf(outfile,"%s_segment%li_calculate_delta_a(_step,cycle%s);\n",indent,segmentNumber,noiseVector);
  fprintf(outfile,"\n");
  
  fprintf(outfile,"%s// c_5==0\n",indent);
  
  fprintf(outfile,"%s// y2=y2+cs_5*a_l\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_check[_i1] += cs[5]*alfield_%s_main[_i1];\n",indent,fieldName,fieldName);

  fprintf(outfile,"%s// a_l= g_1*a_i + g_2*a_k + g_3*a_j + g_4*y_1 +g_5*a_l +g_6*y_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	alfield_%s_main[_i1] = g[1]*aifield_%s_main[_i1] + g[2]*akfield_%s_main[_i1] + g[3]*ajfield_%s_main[_i1] + g[4]*_%s_main[_i1] + g[5]*alfield_%s_main[_i1] + g[6]*_%s_check[_i1] ;\n",indent,fieldName,fieldName,fieldName,fieldName,fieldName,fieldName,fieldName);
  fprintf(outfile,"\n");

  fprintf(outfile,"%s%s += a[6]*_step;\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_l=G[a_l,t+aa_6*dt]\n",indent);
  fprintf(outfile,"%s_segment%li_calculate_delta_a(_step,cycle%s);\n",indent,segmentNumber,noiseVector);
  fprintf(outfile,"\n");
  
  fprintf(outfile,"%s// y1=y1+c_6*a_l\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] += c[6]*alfield_%s_main[_i1];\n",indent,fieldName,fieldName);
  
  fprintf(outfile,"%s// y2=y2+cs_6*a_l\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_check[_i1] += cs[6]*alfield_%s_main[_i1];\n",indent,fieldName,fieldName);
  
  fprintf(outfile,"%s// t->t+dt\n",indent); 
  fprintf(outfile,"%s%s -= a[6]*_step;\n",indent,propDim);
  fprintf(outfile,"\n");
  
  fprintf(outfile,"%s_active_%s_main = _%s_main;\n",indent,fieldName,fieldName);
  fprintf(outfile,"\n"); 
};