/*@@
   @file      sor_confmetric.c
   @date      Tue Sep 26 11:29:18 2000
   @author    Gerd Lanfermann
   @desc 
     Provides sor solver engine routines
   @enddesc 
 @@*/


#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#include "cctk.h"
#include "cctk_Arguments.h"
#include "cctk_Parameters.h"

#include "Boundary.h"
#include "Symmetry.h"
#include "Ell_DBstructure.h"
#include "EllBase.h"

#define SQR(a) ((a)*(a))

int sor_confmetric_3d(cGH *GH, int *MetricPsiI, int conformal,
		       int FieldIndex, int MIndex, int NIndex,
		       CCTK_REAL *AbsTol, CCTK_REAL *RelTol);

 /*@@
   @routine    sor_confmetric
   @date       Tue Sep 26 11:28:08 2000
   @author     Joan Masso, Paul Walker, Gerd Lanfermann
   @desc 
   This is a standalone sor solver engine, 
   called by the wrapper functions
   @enddesc 
   @calls     
   @calledby   
   @history 
 
   @endhistory 

@@*/

int sor_confmetric_3d(cGH *GH, int *MetricPsiI, int conformal,
		       int FieldIndex, int MIndex, int NIndex,
		       CCTK_REAL *AbsTol, CCTK_REAL *RelTol)
{
  DECLARE_CCTK_PARAMETERS  

  int retval = ELL_SUCCESS;
  int timelevel;

  /* The pointer to the metric fields */
  CCTK_REAL *gxx =NULL, *gxy =NULL; 
  CCTK_REAL *gxz =NULL, *gyy =NULL; 
  CCTK_REAL *gyz =NULL, *gzz =NULL; 
  CCTK_REAL *psi =NULL;
  CCTK_REAL *Mlin=NULL, *Nlin=NULL;   
  CCTK_REAL *var =NULL; 

  /* The inverse metric, allocated here */
  CCTK_REAL *uxx=NULL, *uyy=NULL, 
            *uzz=NULL, *uxz=NULL,
            *uxy=NULL, *uyz=NULL;

  /* shortcuts for metric, psi, deltas, etc. */
  CCTK_REAL pm4, p12, det;

  CCTK_REAL dx,dy,dz;
  CCTK_REAL dxdx, dydy, dzdz, 
            dxdy, dxdz, dydz;

  /* Some physical variables */
  int accel_cheb=0, accel_const=0;
  int chebit;
  CCTK_REAL omega, resnorm=0.0, residual=0.0;
  CCTK_REAL glob_residual=0.0, rjacobian=0.0; 
  CCTK_REAL finf;
  CCTK_INT npow;
  CCTK_REAL tol;


  /* Iteration / stepping  variables */
  int sorit; 
  int i,is,ie;
  int j,js,je;
  int k,ks,ke,kstep;
  int nxyz;
  
  /* 19 point stencil index */
  int ijk;
  int ipjk, ijpk, ijkp, imjk, ijmk, ijkm;
  int ipjpk, ipjmk, imjpk, imjmk;
  int ipjkp, ipjkm, imjkp, imjkm;
  int ijpkp, ijpkm, ijmkp, ijmkm;
 
  /* 19 point stencil coefficients */
  CCTK_REAL ac;
  CCTK_REAL ae,aw,an,as,at,ab;
  CCTK_REAL ane, anw, ase, asw, ate, atw, abe, abw;
  CCTK_REAL atn, ats, abn, absol;

  /* Miscellaneous */
  int sum_handle=-1;
  int sw[3], ierr;
  int Mstorage=0, Nstorage=0;
  size_t varsize;
  static int firstcall = 1;
  CCTK_REAL detrec_pm4, sqrtdet;


  /* Get the reduction handle */
  sum_handle = CCTK_ReductionArrayHandle("sum");
  if (sum_handle<0) 
  {
    CCTK_WARN(1,"Cannot get reduction handle for operation >sum<");
  }

  /* IF Robin BCs are set, prepare for a boundary call:
     setup stencil width and get Robin constants (set by the routine
     which is calling the solver interface) */
  if (CCTK_EQUALS(sor_bound,"robin")) 
  { 
    sw[0]=1; 
    sw[1]=1; 
    sw[2]=1;
    
    ierr = Ell_GetRealKey(&finf, "EllLinConfMetric::Bnd::Robin::inf");
    ierr = Ell_GetIntKey (&npow, "EllLinConfMetric::Bnd::Robin::falloff");
  }

  /* Only supports absolute tolerance */
  tol   = AbsTol[0];
  if (CCTK_EQUALS(sor_accel,"const"))
  {
    accel_const = 1;
  }
  else if (CCTK_EQUALS(sor_accel,"cheb"))
  {
    accel_cheb  = 1;
  }

  /* Things to do only once! */
  if (firstcall==1) 
  {
    if (CCTK_Equals(elliptic_verbose, "yes"))
    {
      if (accel_cheb)
      {
        printf("SOR with Chebyshev acceleration with radius of 1\n");
      }
      else if (accel_const)
      {
        printf("SOR with hardcoded omega = 1.8\n");
      }
      else
      {
        printf("SOR with unaccelearted relaxation (omega = 1)\n");
      }
    }
    firstcall = 0;
  }
   

  /* Get the data ptr of these GFs, They all have to be
     on the same timelevel; derive the metric data pointer from 
     the index array. Note the ordering in the metric   */
  timelevel = CCTK_NumTimeLevelsFromVarI (FieldIndex) - 1;
  if (timelevel > 0)
  {
    timelevel--;
  }
  var = (CCTK_REAL*) CCTK_VarDataPtrI(GH, timelevel, FieldIndex);

  timelevel = CCTK_NumTimeLevelsFromVarI (MetricPsiI[0]) - 1;
  if (timelevel > 0)
  {
    timelevel--;
  }
  gxx = (CCTK_REAL*) CCTK_VarDataPtrI(GH, timelevel, MetricPsiI[0]);
  gxy = (CCTK_REAL*) CCTK_VarDataPtrI(GH, timelevel, MetricPsiI[1]);
  gxz = (CCTK_REAL*) CCTK_VarDataPtrI(GH, timelevel, MetricPsiI[2]);
  gyy = (CCTK_REAL*) CCTK_VarDataPtrI(GH, timelevel, MetricPsiI[3]);
  gyz = (CCTK_REAL*) CCTK_VarDataPtrI(GH, timelevel, MetricPsiI[4]);
  gzz = (CCTK_REAL*) CCTK_VarDataPtrI(GH, timelevel, MetricPsiI[5]);

  if (conformal)
  {
    timelevel = CCTK_NumTimeLevelsFromVarI (MetricPsiI[6]) - 1;
    if (timelevel > 0)
    {
      timelevel--;
    }
    psi = (CCTK_REAL*) CCTK_VarDataPtrI(GH, timelevel, MetricPsiI[6]);
  }

  /* if we have a negative index for M/N, this GF is not needed, 
     there for don't even look for it. when index positive,
     set flag Mstorage=1, dito for N */
  if (MIndex>=0)  
  { 
    timelevel = CCTK_NumTimeLevelsFromVarI (MIndex) - 1;
    if (timelevel > 0)
    {
      timelevel--;
    }
    Mlin = (CCTK_REAL*) CCTK_VarDataPtrI(GH,timelevel,MIndex);
    Mstorage = 1;
  }
  if (NIndex>=0) 
  {
    timelevel = CCTK_NumTimeLevelsFromVarI (NIndex) - 1;
    if (timelevel > 0)
    {
      timelevel--;
    }
    Nlin = (CCTK_REAL*) CCTK_VarDataPtrI(GH,timelevel,NIndex);
    Nstorage = 1;
  }

  /* Shortcuts */
  dx   = GH->cctk_delta_space[0];
  dy   = GH->cctk_delta_space[1];
  dz   = GH->cctk_delta_space[2];
  nxyz = GH->cctk_lsh[0]*GH->cctk_lsh[1]*GH->cctk_lsh[2];


  /* Allocate the inverse metric */
  varsize = (size_t)CCTK_VarTypeSize(CCTK_VarTypeI(FieldIndex))*nxyz;
  if (!(uxx = (CCTK_REAL*) malloc(varsize))) CCTK_WARN(0,"Allocation failed");
  if (!(uyy = (CCTK_REAL*) malloc(varsize))) CCTK_WARN(0,"Allocation failed");
  if (!(uzz = (CCTK_REAL*) malloc(varsize))) CCTK_WARN(0,"Allocation failed");
  if (!(uxy = (CCTK_REAL*) malloc(varsize))) CCTK_WARN(0,"Allocation failed");
  if (!(uxz = (CCTK_REAL*) malloc(varsize))) CCTK_WARN(0,"Allocation failed");
  if (!(uyz = (CCTK_REAL*) malloc(varsize))) CCTK_WARN(0,"Allocation failed");

  /* PreCalc the differential coeff. */
  dxdx = 1.0/(2.0*dx*dx);
  dydy = 1.0/(2.0*dy*dy);
  dzdz = 1.0/(2.0*dz*dz);
  dxdy = 1.0/(4.0*dx*dy);
  dxdz = 1.0/(4.0*dx*dz);
  dydz = 1.0/(4.0*dy*dz);

  /* Calculate the inverse metric */

  for (ijk=0;ijk<nxyz;ijk++) 
  {

    /* determinant */
    det = -(SQR(gxz[ijk])*gyy[ijk]) + 
      2*gxy[ijk]*gxz[ijk]*gyz[ijk] - 
      gxx[ijk]*SQR(gyz[ijk])  -
      SQR(gxy[ijk])*gzz[ijk] + 
      gxx[ijk]*gyy[ijk]*gzz[ijk];
    
    if (conformal) 
    {
      pm4 = 1.0/pow(psi[ijk],4.0);
      p12 = pow(psi[ijk],12.0);
    } 
    else 
    {
      pm4 = 1.0;
      p12 = 1.0;
    }


    /* try to avoid constly division */
    detrec_pm4 = 1.0/det*pm4;
    sqrtdet    = sqrt(det);

    uxx[ijk]=(-SQR(gyz[ijk])     + gyy[ijk]*gzz[ijk])*detrec_pm4;
    uxy[ijk]=( gxz[ijk]*gyz[ijk] - gxy[ijk]*gzz[ijk])*detrec_pm4;
    uxz[ijk]=(-gxz[ijk]*gyy[ijk] + gxy[ijk]*gyz[ijk])*detrec_pm4;
    uyy[ijk]=(-SQR(gxz[ijk])     + gxx[ijk]*gzz[ijk])*detrec_pm4;
    uyz[ijk]=( gxy[ijk]*gxz[ijk] - gxx[ijk]*gyz[ijk])*detrec_pm4;
    uzz[ijk]=(-SQR(gxy[ijk])     + gxx[ijk]*gyy[ijk])*detrec_pm4;
 	
    det    = det*p12;
    sqrtdet= sqrt(det);
   
    /* Rescaling for the uppermetric solver */
    if (Mstorage) Mlin[ijk] = Mlin[ijk]*sqrt(det);
    if (Nstorage) Nlin[ijk] = Nlin[ijk]*sqrt(det);
    
    uxx[ijk]=uxx[ijk]*dxdx*sqrtdet;
    uyy[ijk]=uyy[ijk]*dydy*sqrtdet;
    uzz[ijk]=uzz[ijk]*dzdz*sqrtdet;
    uxy[ijk]=uxy[ijk]*dxdy*sqrtdet;
    uxz[ijk]=uxz[ijk]*dxdz*sqrtdet;
    uyz[ijk]=uyz[ijk]*dydz*sqrtdet;

  }	

  /*$for (k=1;k<GH->cctk_lsh[2]-1;k++) {
    for (j=1;j<GH->cctk_lsh[1]-1;j++)   
    {
      for (i=1;i<GH->cctk_lsh[0]-1;i++)   
      {
	ijk   = CCTK_GFINDEX3D(GH,i,j,k); 
	printf("U G : %d %d %d   %7.8g  %7.8g  %7.8g \n",
	       i,j,k,uxx[ijk],uyy[ijk],gxx[ijk]);
      }
    }
  }$*/
	

  /* iteration boundaries (ks set inside loop)*/
  is = 1;
  js = 1;
  ie = GH->cctk_lsh[0]-1;
  je = GH->cctk_lsh[1]-1;
  ke = GH->cctk_lsh[2]-1;
  kstep = 2;

  /* start at 1 for historic (Fortran) reasons */
  for (sorit=1; sorit<=maxit; sorit++) {
    
    omega     =  1.0;
    rjacobian =  1.0;

    if (accel_cheb) 
    {
      for (chebit=2;chebit<sorit;chebit++)
      {
	omega = 1.0/(1.0 - 0.25*rjacobian*rjacobian*omega);
      }
    }
    if (accel_const)
    {
      omega = 1.8;
    }

    resnorm = 0.0;

    ks = (sorit%2)+1;
    if (GH->cctk_lsh[2]==3)
      ks = 1;

    for (k=ks;k<ke;k+=kstep) 
    {
      for (j=js;j<je;j++)     
      {
	for (i=is;i<ie;i++)    
        {

	  ijk   = CCTK_GFINDEX3D(GH,i  ,j  ,k  );
	  
	  ipjk  = CCTK_GFINDEX3D(GH,i+1,j  ,k  );
	  imjk  = CCTK_GFINDEX3D(GH,i-1,j  ,k  );
	  ijpk  = CCTK_GFINDEX3D(GH,i  ,j+1,k  );
	  ijmk  = CCTK_GFINDEX3D(GH,i  ,j-1,k  );
	  ijkp  = CCTK_GFINDEX3D(GH,i  ,j  ,k+1);
	  ijkm  = CCTK_GFINDEX3D(GH,i  ,j  ,k-1);
	  
	  ipjpk = CCTK_GFINDEX3D(GH,i+1,j+1,k  );
	  ipjmk = CCTK_GFINDEX3D(GH,i+1,j-1,k  );
	  imjpk = CCTK_GFINDEX3D(GH,i-1,j+1,k  );
	  imjmk = CCTK_GFINDEX3D(GH,i-1,j-1,k  );
	  
	  ijpkp = CCTK_GFINDEX3D(GH,i  ,j+1,k+1);
	  ijpkm = CCTK_GFINDEX3D(GH,i  ,j+1,k-1);
	  ijmkp = CCTK_GFINDEX3D(GH,i  ,j-1,k+1);
	  ijmkm = CCTK_GFINDEX3D(GH,i  ,j-1,k-1);
	  
	  ipjkp = CCTK_GFINDEX3D(GH,i+1,j  ,k+1);
	  ipjkm = CCTK_GFINDEX3D(GH,i+1,j  ,k-1);
	  imjkp = CCTK_GFINDEX3D(GH,i-1,j  ,k+1);
	  imjkm = CCTK_GFINDEX3D(GH,i-1,j  ,k-1);
	  

	  ac = -1.0*uxx[ipjk] - 2.0*uxx[ijk] - 1.0*uxx[imjk]
	    -1.0*uyy[ijpk] - 2.0*uyy[ijk] - 1.0*uyy[ijmk]
	    -1.0*uzz[ijkp] - 2.0*uzz[ijk] - 1.0*uzz[ijkm];
	  
	  
	  if (Mstorage) 
	  {
	    ac += Mlin[ijk];
	  }

	  ae  = uxx[ipjk]+uxx[ijk];
	  aw  = uxx[imjk]+uxx[ijk];
	  an  = uyy[ijpk]+uyy[ijk];
	  as  = uyy[ijmk]+uyy[ijk];
	  at  = uzz[ijkp]+uzz[ijk];
	  ab  = uzz[ijkm]+uzz[ijk];
	  
	  ane = uxy[ijpk]+uxy[ipjk];
	  anw =-uxy[imjk]-uxy[ijpk];
	  ase =-uxy[ipjk]-uxy[ijmk];
	  asw = uxy[imjk]+uxy[ijmk];
	  
	  ate = uxz[ijkp]+uxz[ipjk];
	  atw =-uxz[imjk]-uxz[ijkp];
	  abe =-uxz[ipjk]-uxz[ijkm];
	  abw = uxz[imjk]+uxz[ijkm];
	  
	  atn = uyz[ijpk]+uyz[ijkp];               
	  ats =-uyz[ijkp]-uyz[ijmk];
	  abn =-uyz[ijkm]-uyz[ijpk];
	  absol = uyz[ijkm]+uyz[ijmk];
	  
	  residual = ac * var[ijk]
	    + ae *var[ipjk]  +  aw*var[imjk]
	    + an *var[ijpk]  +  as*var[ijmk]
	    + at *var[ijkp]  +  ab*var[ijkm];
	  
	  residual = residual 
	    + ane*var[ipjpk] + anw*var[imjpk]; 
	  
	  residual = residual 
	    + ase*var[ipjmk] + asw*var[imjmk];
	  
	  residual = residual    
	    + ate*var[ipjkp] + atw*var[imjkp] 
	    + abe*var[ipjkm] + abw*var[imjkm]
	    + atn*var[ijpkp] + ats*var[ijmkp]
	    + abn*var[ijpkm] + absol*var[ijmkm];
	  
	  if (Nstorage)
	  {
	    residual +=Nlin[ijk];
	  }
    
	  resnorm  = resnorm + fabs(residual);
	  var[ijk] = var[ijk] - omega*residual/ac; 
	}
      }
    }

    /* reduction operation on processor-local residual values */
    ierr = CCTK_ReduceLocScalar(GH, -1, sum_handle,
				&resnorm, &glob_residual, CCTK_VARIABLE_REAL); 
    if (ierr<0) 
    {
      CCTK_WARN(1,"Reduction of residual  failed");
    }

    glob_residual = glob_residual /  
      (GH->cctk_gsh[0]*GH->cctk_gsh[1]*GH->cctk_gsh[2]);

    /* apply symmetry boundary conditions within loop */    
    if (CartSymVI(GH,FieldIndex)<0) 
    {
      CCTK_WARN(1,"CartSymVI failed in EllSOR loop");
    }

    /* apply boundary conditions within loop */
    if (CCTK_EQUALS(sor_bound,"robin"))
    {
      ierr = BndRobinVI(GH, sw, finf, npow,  FieldIndex);
    }

    /* synchronization of grid variable */
    CCTK_SyncGroupWithVarI(GH, FieldIndex);

    /* Leave iteration loop if tolerance criterium is met */
    if (glob_residual<tol)
    {
      break;
    }
  }

  /* Information for the user if the solve did not converge within 
     the given constraints of max.iteration and tolerance */
  if (residual>tol) 
  {
    retval = ELL_NOCONVERGENCE;
    CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING,
	       "SOR SOLVER DID NOT CONVERGE within given "
	       "tolerance/max.number of iterations.\n "
	       "max. iteration %d   residual (tolerance): %g (%g)\n",
	       maxit, glob_residual, tol );
  }

  if (uxx) free(uxx); 
  if (uyy) free(uyy); 
  if (uzz) free(uzz);
  if (uxy) free(uxy); 
  if (uxz) free(uxz); 
  if (uyz) free(uyz);

  return retval;
}