truncation.C

//------------------------------------------------------------------------
// Routine to carry out a truncation error estimate 
// in first order form
//------------------------------------------------------------------------

#include "DAGH.h"
#include "tportamr.h"
#include "tportfortran.h"

#define VizServer	0	// Who will write to stdout

void wave_trunc_est(GridHierarchy &GH,
      INTEGER const t,
      INTEGER const l,
      GridFunction(DIM)<GFTYPE> &U,
      GridFunction(DIM)<GFTYPE> &trunc_err)
{ 

      GFTYPE normu, inormu, tnormu, enorm;

      GridFunction(DIM)<GFTYPE> truncU("truncU",U,t,l, 
                                       DAGHComm, 
                                       DAGHHasShadow,
                                       DAGHNoBoundary,
				       DAGHNoAdaptBoundary);

//------------------------------------------------------------------------
// Subtract the DAGH_Main grid from shadow pointwise
//------------------------------------------------------------------------
      forall(truncU, t, l, c)
         truncU(t,l,c,DAGH_Shadow) = U(t,l,c,DAGH_Shadow);
         truncU(t,l,c,DAGH_Shadow) -= U(t,l,c,DAGH_Main);
      end_forall

//------------------------------------------------------------------------
// "Relativize" wrt to norms of U, (DAGH_Main)
//
// 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.
//------------------------------------------------------------------------
      normu = Norm2(U, t, l, DAGH_Main);
      if (normu == 0.0) {
         inormu = 1.0;
      } else {
         inormu = 1.0 / (normu*normu);
      }
      
//------------------------------------------------------------------------
// Take a Pythogorean norm of the truncation error
// trunc_err must have the same bbox properties as truncU.
//------------------------------------------------------------------------
      forall(trunc_err, t, l, c)
        f_pythnorm(FBA(trunc_err(t,l,c,DAGH_Shadow)), 
                   FDA(truncU(t,l,c,DAGH_Shadow)), 
                   FDA(trunc_err(t,l,c,DAGH_Shadow)),
                   &inormu);
      end_forall

      Sync(trunc_err, t, l, DAGH_Shadow);

//------------------------------------------------------------------------
// Prolong DAGH_Shadow trunc_err to DAGH_Main trunc_err
//------------------------------------------------------------------------
      INTEGER myargc = 0; GFTYPE myargs[1];
      forall(trunc_err,t,l,c)
        f_prolong_amr(FDA(trunc_err(t,l,c,DAGH_Shadow)),
                      FBA(trunc_err(t,l,c,DAGH_Shadow)),
                      FDA(trunc_err(t,l,c,DAGH_Main)), 
                      FBA(trunc_err(t,l,c,DAGH_Main)), 
                      FBA(trunc_err(t,l,c,DAGH_Main)), 
                      myargs, &myargc); 
      end_forall

      Sync(trunc_err, t, l, DAGH_Main);
} /* end wave_trunc_est */