/*------------------------------------------------------------------------*/
// truncation.C
//
// Truncation error estimators
/*------------------------------------------------------------------------*/
#include "DAGH.h"
#include "waveamr2d.h"
#include "wavefortran2d.h"

#define VizServer 0

// Routine to compute truncation error estimate for single grid function
// truncgf has to be declared with the same size as shadowgf...
// that is declare as 
// GridFunction(DIM)<DOUBLE> truncgf("truncgf",shadowgf);
//

void truncation_est(GridHierarchy &GH,
      INTEGER const t,
      INTEGER const l,
      GridFunction(DIM)<DOUBLE> &maingf,
      GridFunction(DIM)<DOUBLE> &shadowgf,
      GridFunction(DIM)<DOUBLE> &truncgf)
{ /* truncation_est */
   INTEGER t_sten = 0, s_sten = 1;

   forall(truncgf, t, l, c)
      truncgf(t,l,c) -= maingf(t,l,c);
   end_forall

} /* end truncation_est */

/*-----------------------------------------------------------------------*/
// Routine to carry out a truncation error estimate for the 2D wave equation
// in first order form
//
// What has to be passed in are A, B, C and phi
// 
void wave_trunc_est(GridHierarchy &GH,
      INTEGER const t,
      INTEGER const l,
      GridFunction(DIM)<DOUBLE> &A,
      GridFunction(DIM)<DOUBLE> &B,
      GridFunction(DIM)<DOUBLE> &C,
      GridFunction(DIM)<DOUBLE> &trunc_err)
{ /* wave_trunc_est */

      DOUBLE norma, normb, normc, enorm;
      DOUBLE inorma, inormb, inormc;
      DOUBLE tnorma, tnormb, tnormc;


// Declare grid functions for storing truncation errors for grid functions
// A, B, and C. Ensure that they share the same properties as the SHADOW's.
      GridFunction(DIM)<DOUBLE> truncA("trunca",A,t,l, COMM, HAS_SHADOW);
      GridFunction(DIM)<DOUBLE> truncB("truncb",B,t,l, COMM, HAS_SHADOW);
      GridFunction(DIM)<DOUBLE> truncC("truncc",C,t,l, COMM, HAS_SHADOW);


// Now subtract the MAIN grid off pointwise

      forall(truncA, t, l, c)
         truncA(t,l,c,SHADOW) = A(t,l,c,SHADOW);
      end_forall

      forall(truncB, t, l, c)
         truncB(t,l,c,SHADOW) = B(t,l,c,SHADOW);
      end_forall

      forall(truncC, t, l, c)
         truncC(t,l,c,SHADOW) = C(t,l,c,SHADOW);
      end_forall

      forall(truncA, t, l, c)
         truncA(t,l,c,SHADOW) -= A(t,l,c,MAIN);
      end_forall

      forall(truncB, t, l, c)
         truncB(t,l,c,SHADOW) -= B(t,l,c,MAIN);
      end_forall

      forall(truncC, t, l, c)
         truncC(t,l,c,SHADOW) -= C(t,l,c,MAIN);
      end_forall


// "Relativize" wrt to norms of A, B, C (MAIN)

      norma = Norm2(A, t, l, MAIN);
      normb = Norm2(B, t, l, MAIN);
      normc = Norm2(C, t, l, MAIN);

      DOUBLE shnorma = Norm2(A, t, l, SHADOW);
      DOUBLE shnormb = Norm2(B, t, l, SHADOW);
      DOUBLE shnormc = Norm2(C, t, l, SHADOW);

      tnorma = Norm2(truncA, t, l, SHADOW);
      tnormb = Norm2(truncB, t, l, SHADOW);
      tnormc = Norm2(truncC, t, l, SHADOW);

      if(MY_PROC(GH) == 0){
        cout << "Truncation error: " 
             << tnorma << " " << tnormb << " " << tnormc << endl << flush ;
        cout << "Norm of A,B,C " 
             << norma << " " << normb << " " << normc << endl << flush ;
        cout << "Shadow Norm of A,B,C " 
             << shnorma << " " << shnormb << " " << shnormc << endl << flush ;
      }

//   Check if the norms are too small... If they are then we assume that
// the grid functions consist of small quantities and do not divide by
// the norm.
      
      if(norma == 0.0){
         inorma = 1.0;
      }else{
         inorma = 1.0 / (norma*norma);
      }
      
      if(normb == 0.0){
         inormb = 1.0;
      }else{
         inormb = 1.0 / (normb*normb);
      }

      if(normc == 0.0){
         inormc = 1.0;
      }else{
         inormc = 1.0 / (normc*normc);
      }

      if(MY_PROC(GH) == 0){
         cout << inorma << " " << inormb << " " << inormc << endl<< flush;
      }

// Take a Pythogorean norm of the three truncation errors 
// trunc_err must have the same bbox properties as truncA, truncB and truncC.

      forall(trunc_err, t, l, c)
        f_pythnorm(FBA(trunc_err(t,l,c,SHADOW)), 
                   FDA(truncA(t,l,c,SHADOW)),FDA(truncB(t,l,c,SHADOW)), 
                   FDA(truncC(t,l,c,SHADOW)),FDA(trunc_err(t,l,c,SHADOW)),
                   &inorma,&inormb,&inormc);
      end_forall


      Sync(trunc_err, t, l, SHADOW);

      enorm = Norm2(trunc_err, t, l, SHADOW);

      // Prolong SHADOW trunc_err to MAIN trunc_err
      forall(trunc_err,t,l,c)
        f_prolong_amr(FDA(trunc_err(t,l,c,SHADOW)),FBA(trunc_err(t,l,c,SHADOW)),
                      FDA(trunc_err(t,l,c,MAIN)), FBA(trunc_err(t,l,c,MAIN)), 
                      FBA(trunc_err(t,l,c,MAIN)),0, 0); 
      end_forall

      Sync(trunc_err, t, l, MAIN);



      cout << "L2 Norm of truncation error (SHADOW) at Level "  << l << " " 
           << enorm << endl<< flush ;
      cout << "L2 Norm of truncation error (MAIN) at Level "  << l << " " 
           << Norm2(trunc_err, t, l, MAIN) << endl<< flush ;

      cout << "MaxVal of truncation error (MAIN) at Level "  << l << " "  
           << MaxVal(trunc_err, t, l, MAIN) << endl<< flush ;
      cout << "MaxVal of truncation error (SHADOW) at Level "  << l << " "
           << MaxVal(trunc_err, t, l, SHADOW) << endl<< flush ;

      // IO for one processor
      if (1) {
      char fname[40];
      static int cnt = 0;
      sprintf(fname, "truncLevel%derr.dat%d",l,cnt);
      cnt++;
      writeslice(MAIN, fname, GH, trunc_err, t, l, VizServer);
      }
      
} /* end wave_trunc_est */