/*************************************************************************/
/* Copyright (C) 1995-1997 Marlis Hochbruck, Christian Lubich, Hubert    */
/* Selhofer. All rights reserved.                                        */
/*                                                                       */
/* C code written by Mathias Froehlich.                                  */
/*                                                                       */
/* This code is part of a copyrighted package. For details, see the file */
/* "COPYING" in the top-level directory.                                 */
/*************************************************************************/

/************************************************************************/
/* The brusselator example in two dimensions.                           */
/* The N defined below is the number of grid points per dimension.      */
/* The resulting problem dimension is 2*N*N.                            */
/* For details about the equation see                                   */
/* Hairer, Noersett, Wanner: Solving Ordinary Differential              */
/* Equations II, p 159. Springer-Verlag 1991.                           */
/************************************************************************/

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "ftypes.h"  /* needed for the types fdouble,finteger ...    */
#include "fblas.h"   /* all the F77 blas - routines are defined here */
#include "exp4.h"    /* Exp4 User Routines ...                       */

#define N         20
#define ALPHA     0.02
/* #define N         70 */
/* #define ALPHA     0.2 */

#define INITIAL_VALUES Peaks
/* #define INITIAL_VALUES Planes */

#define   Nm1   (N-1)
#define   Nm2   (N-2)
#define   Nm3   (N-3)
#define   NN    (N*N)

#define   DxDx1Alpha (ALPHA*(Nm1)*(Nm1))
#define   DxDx2Alpha (2.0*ALPHA*(Nm1)*(Nm1))
#define   DxDx4Alpha (4.0*ALPHA*(Nm1)*(Nm1))

/* static data for the righthandside function */
static fdouble *J11;
static fdouble *J12;
static fdouble *J21;
static fdouble *J22;

/* the selection of the initial conditions */
typedef enum { Peaks, Planes } Bruss2dStart;

/* the righthandside function itself, calculates fy:=f(y) */
void Bruss2d(finteger neq,fdouble *fy,fdouble *y)
{
#define U(I,J)  y[(I)+(J)*N]
#define V(I,J)  y[NN+(I)+(J)*N]
#define dU(I,J) fy[(I)+(J)*N]
#define dV(I,J) fy[NN+(I)+(J)*N]
  finteger i,j;

  for (j=0;j<N;j++)
    for (i=0;i<N;i++)
    {
      fdouble u = U(i,j);
      fdouble v = V(i,j);
      fdouble uuv = u*u*v;
      fdouble du = 1.0 + uuv - 4.4*u  - DxDx4Alpha*u;
      fdouble dv = 3.4*u - uuv  - DxDx4Alpha*v;

      if (i == 0) {
	du += DxDx2Alpha*U(1,j);
	dv += DxDx2Alpha*V(1,j);
      } else if (i == Nm1) {
	du += DxDx2Alpha*U(Nm2,j);
	dv += DxDx2Alpha*V(Nm2,j);
      } else {
	du += DxDx1Alpha*(U(i-1,j)+U(i+1,j));
	dv += DxDx1Alpha*(V(i-1,j)+V(i+1,j));
      }

      if (j == 0) {
	du += DxDx2Alpha*U(i,1);
	dv += DxDx2Alpha*V(i,1);
      } else if (j == Nm1) {
	du += DxDx2Alpha*U(i,Nm2);
	dv += DxDx2Alpha*V(i,Nm2);
      } else {
	du += DxDx1Alpha*(U(i,j-1)+U(i,j+1));
	dv += DxDx1Alpha*(V(i,j-1)+V(i,j+1));
      }

      dU(i,j) = du;
      dV(i,j) = dv;
    }
#undef U
#undef V
#undef dU
#undef dV
}

/* the Jacobian function, calculates the y dependant values of the Jacobian */
void Bruss2dNewJac(finteger neq,fdouble *y)
{
#define U(I,J)     y[(I)+(J)*N]
#define V(I,J)     y[NN+(I)+(J)*N]
  finteger i,j;

  for (j=0;j<N;j++)
    for (i=0;i<N;i++)
    {
      fdouble u = U(i,j);
      fdouble v = V(i,j);
      fdouble uv = 2.0*u*v;
      fdouble uu = u*u;
      J11[i+N*j] = uv - 4.4  - DxDx4Alpha;
      J12[i+N*j] = uu;
      J22[i+N*j] = - uu - DxDx4Alpha;
      J21[i+N*j] = 3.4 - uv;      
    }
#undef U
#undef V
}

/* the multiplication routine, res := Jacobian*vect */
void Bruss2dJacMultV(finteger neq,fdouble *res,fdouble *dummy,fdouble* vect)
{
#define U(I,J)     vect[(I)+(J)*N]
#define V(I,J)     vect[NN+(I)+(J)*N]
#define URes(I,J)  res[(I)+(J)*N]
#define VRes(I,J)  res[NN+(I)+(J)*N]
  finteger i,j;

  for (j=0;j<N;j++)
    for (i=0;i<N;i++)
    {
      fdouble u = U(i,j);
      fdouble v = V(i,j);
      fdouble du = J11[i+N*j]*u + J12[i+N*j]*v;
      fdouble dv = J21[i+N*j]*u + J22[i+N*j]*v;

      if (i == 0) {
	du += DxDx2Alpha*U(1,j);
	dv += DxDx2Alpha*V(1,j);
      } else if (i == Nm1) {
	du += DxDx2Alpha*U(Nm2,j);
	dv += DxDx2Alpha*V(Nm2,j);
      } else {
	du += DxDx1Alpha*(U(i-1,j)+U(i+1,j));
	dv += DxDx1Alpha*(V(i-1,j)+V(i+1,j));
      }

      if (j == 0) {
	du += DxDx2Alpha*U(i,1);
	dv += DxDx2Alpha*V(i,1);
      } else if (j == Nm1) {
	du += DxDx2Alpha*U(i,Nm2);
	dv += DxDx2Alpha*V(i,Nm2);
      } else {
	du += DxDx1Alpha*(U(i,j-1)+U(i,j+1));
	dv += DxDx1Alpha*(V(i,j-1)+V(i,j+1));
      }

      URes(i,j) = du;
      VRes(i,j) = dv;
    }
#undef U
#undef V
#undef URes
#undef VRes
}

/* the transpose multiplication routine, res := Jacobian'*vect */
void Bruss2dJacTMultV(finteger neq,fdouble *res,fdouble *dummy,fdouble* vect)
{
#define U(I,J)     vect[(I)+(J)*N]
#define V(I,J)     vect[NN+(I)+(J)*N]
#define URes(I,J)  res[(I)+(J)*N]
#define VRes(I,J)  res[NN+(I)+(J)*N]
  finteger i,j;

  for (j=0;j<N;j++)
    for (i=0;i<N;i++)
    {
      fdouble u = U(i,j);
      fdouble v = V(i,j);
      fdouble du = J11[i+N*j]*u + J21[i+N*j]*v;
      fdouble dv = J12[i+N*j]*u + J22[i+N*j]*v;

      if (i == 0) {
	du += DxDx1Alpha*U(1,j);
	dv += DxDx1Alpha*V(1,j);
      } else if (i == 1) {
	du += DxDx1Alpha*(2.0*U(0,j)+U(2,j));
	dv += DxDx1Alpha*(2.0*V(0,j)+V(2,j));
      } else if (i == Nm2) {
	du += DxDx1Alpha*(U(Nm3,j)+2.0*U(Nm1,j));
	dv += DxDx1Alpha*(V(Nm3,j)+2.0*V(Nm1,j));
      } else if (i == Nm1) {
	du += DxDx1Alpha*U(Nm2,j);
	dv += DxDx1Alpha*V(Nm2,j);
      } else {
	du += DxDx1Alpha*(U(i-1,j)+U(i+1,j));
	dv += DxDx1Alpha*(V(i-1,j)+V(i+1,j));
      }

      if (j == 0) {
	du += DxDx1Alpha*U(i,1);
	dv += DxDx1Alpha*V(i,1);
      } else if (j == 1) {
	du += DxDx1Alpha*(2.0*U(i,0)+U(i,2));
	dv += DxDx1Alpha*(2.0*V(i,0)+V(i,2));
      } else if (j == Nm2) {
	du += DxDx1Alpha*(U(i,Nm3)+2.0*U(i,Nm1));
	dv += DxDx1Alpha*(V(i,Nm3)+2.0*V(i,Nm1));
      } else if (j == Nm1) {
	du += DxDx1Alpha*U(i,Nm2);
	dv += DxDx1Alpha*V(i,Nm2);
      } else {
	du += DxDx1Alpha*(U(i,j-1)+U(i,j+1));
	dv += DxDx1Alpha*(V(i,j-1)+V(i,j+1));
      }

      URes(i,j) = du;
      VRes(i,j) = dv;
    }
#undef U
#undef V
#undef URes
#undef VRes
}

/* initialisation function, allocates memory needed and sets the
 * initial conditions and the parameters of the problem.
 * (alpha: diffusion coeficient, n: # of gridpoints in each direction) */
fdouble *InitBruss2d()
{
  fdouble *y0;

  if (NULL == (y0 = malloc(6*NN*sizeof(fdouble)))) {
    fprintf(stderr,"InitBruss2d: can't get memory for working!\n");
    exit(0);
  }
  J11 = y0+2*NN;
  J12 = J11+NN;
  J22 = J12+NN;
  J21 = J22+NN;

  /* generating start values y0 */
  switch (INITIAL_VALUES) {
  case Peaks:
    {
#define X(I) (-3.0+6.0/((fdouble)Nm1)*I)
#define Y(I) (-3.0+6.0/((fdouble)Nm1)*I)
      finteger i,j;
      for (i=0;i<N;i++)
	for (j=0;j<N;j++)
	{
	  y0[i*N+j]=3.0*(1.0-X(i))*(1.0-X(i))*
	    exp(-X(i)*X(i)-(Y(j)+1.0)*(Y(j)+1.0));
	  y0[i*N+j]-= 10.0*(X(i)/5.0 - X(i)*X(i)*X(i) - 
		Y(j)*Y(j)*Y(j)*Y(j)*Y(j))*exp(-X(i)*X(i)-Y(j)*Y(j));
	  y0[i*N+j]-=1.0/3.0*exp(-(X(i)+1.0)*(X(i)+1.0)- Y(j)*Y(j));
	}
      for (i=NN;i<2*NN;i++)
	y0[i]=0.0;

      break;
#undef X
#undef Y
    }
  case Planes:
    {
      finteger i,j;

      for (i=0;i<N;i++)
	for (j=0;j<N;j++)
	{
	  y0[i*N+j]    = 0.5 + 1.0/Nm1*j;
	  y0[NN+i*N+j] = 1.0 + 5.0/Nm1*i;
	}

      break;
    }
  }

  return y0;
}

int main()
{
  /* the start time of integration */
  fdouble t=0.0;
  /* space for the initial data and the returning solution vector */
  fdouble *y0;
  /* number of equations that is the dimension of im(f) */
  finteger neq = 2*N*N;

  /* initialize the right-hand side function and the jacobian function */
  /* generating start values y0 */
  y0 = InitBruss2d();

  printf("\n#################################################################"
	"\nExp4 Example: The Brusselator example, %d - dimensional Problem\n"
	"#################################################################\n"
	"\nSimple integrator call!\n"
	"\nPress return please ...\n",(int)neq);
  getchar();

  /* Documentation of the functions below in ../exp4/exp4.h */

  /* set Atol */
/*   DExp4SetAtol(Scalar,1e-8); */

  /* set Rtol */
/*   DExp4SetRtol(Scalar,1e-6); */

  /* call the integrator */
  DExp4(&t,1e1,y0,neq,&Bruss2d,&Bruss2dNewJac,&Bruss2dJacMultV,Exp4False);

  /* free memory */
  free(y0);

  return 1;
}