
#include <new>
#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <string>
using std::ios;

#include <unistd.h>

#include <AMReX_REAL.H>
#include <AMReX_Box.H>
#include <AMReX_FArrayBox.H>
#include <AMReX_ParmParse.H>
#include <AMReX_ParallelDescriptor.H>
#include <AMReX_DataServices.H>
#include <AMReX_Utility.H>
#include <AMReX_VisMF.H>


#define GARBAGE 666.e+40
using namespace amrex;
static
void
PrintUsage (const char* progName)
{
    std::cout << '\n';
    std::cout << "Usage:" << '\n';
    std::cout << progName << '\n';
    std::cout << "    infile = inputFileName" << '\n';
    std::cout << "      norm = integer norm (Ie. default is 2 for L2 norm)" << '\n';
    std::cout << "   [-help]" << '\n';
    std::cout << "   [-verbose]" << '\n';
    std::cout << '\n';
    exit(1);
}

int
main (int   argc,
      char* argv[])
{
    amrex::Initialize(argc,argv);
    {

        if (argc == 1)
            PrintUsage(argv[0]);

        ParmParse pp;

        const std::string farg = amrex::get_command_argument(1);
        if (farg == "-h" || farg == "--help")
            PrintUsage(argv[0]);

        FArrayBox::setFormat(FABio::FAB_IEEE_32);
        //
        // Scan the arguments.
        //
        std::string iFile;

        bool verbose = false;
        if (pp.contains("verbose"))
        {
            verbose = true;
            AmrData::SetVerbose(true);
        }
        pp.query("infile", iFile);
        if (iFile.empty())
            amrex::Abort("You must specify `infile'");

        int norm = 2;
        pp.query("norm", norm);

        DataServices::SetBatchMode();
        Amrvis::FileType fileType(Amrvis::NEWPLT);

        DataServices dataServicesC(iFile, fileType);

        if (!dataServicesC.AmrDataOk())
            amrex::Abort("ERROR: Dataservices not OK");


        //
        // Generate AmrData Objects
        //
        AmrData& amrDataI = dataServicesC.AmrDataRef();

        //
        // Initial Tests
        //
        int nComp       = amrDataI.NComp();
        int finestLevel = amrDataI.FinestLevel();
        const Vector<std::string>& derives = amrDataI.PlotVarNames();
        Vector<int> destComps(nComp);
        for (int i = 0; i < nComp; i++)
            destComps[i] = i;


        //
        // Compute the error
        //
        Vector<MultiFab*> error(finestLevel+1);

        if (ParallelDescriptor::IOProcessor())
            std::cout << "Level  L"<< norm << " norm of Error in Each Component" << std::endl
                      << "-----------------------------------------------" << std::endl;

        for (int iLevel = 0; iLevel <= finestLevel; ++iLevel)
        {
            const BoxArray& baI = amrDataI.boxArray(iLevel);

            DistributionMapping dm(baI);
            MultiFab dataI(baI, dm, nComp, 0);

            amrDataI.FillVar(dataI, iLevel, derives, destComps);

            //
            // Output Statistics
            //
            if (ParallelDescriptor::IOProcessor())
                std::cout << "  " << iLevel << "    ";

            Vector<Real> norms(nComp);
            for (int iComp = 0; iComp < nComp; iComp++)
                norms[iComp] = 0.0;

            for (MFIter mfi(dataI); mfi.isValid(); ++mfi)
            {
                for (int iComp = 0; iComp < nComp; iComp++)
                {
                    Real grdL2 = dataI[mfi].norm(norm, iComp, 1);

                    if (norm != 0)
                    {
                        norms[iComp] = norms[iComp] + pow(grdL2, norm);
                    }
                    else
                    {
                        norms[iComp] = std::max(norms[iComp], grdL2);
                    }

                }
            }


#ifdef BL_USE_MPI
            MPI_Datatype datatype = mpi_data_type(norms.dataPtr());
            if (ParallelDescriptor::IOProcessor())
            {
                Vector<Real> tmp(nComp);
                for (int proc = 0; proc < ParallelDescriptor::NProcs(); proc++)
                    if (proc != ParallelDescriptor::IOProcessorNumber())
                    {
                        MPI_Status stat;
                        int rc = MPI_Recv(tmp.dataPtr(), nComp, datatype,
                                          MPI_ANY_SOURCE, proc, ParallelDescriptor::Communicator(),
                                          &stat);

                        if (rc != MPI_SUCCESS)
                            ParallelDescriptor::Abort(rc);

                        for (int iComp = 0; iComp < nComp; iComp++)
                            if (norm != 0)
                            {
                                norms[iComp] = norms[iComp] + tmp[iComp];
                            }
                            else
                            {
                                norms[iComp] = Max(norms[iComp], tmp[iComp]);
                            }
                    }
            }
            else
            {
                int rc = MPI_Send(norms.dataPtr(), nComp, datatype,
                                  ParallelDescriptor::IOProcessorNumber(),
                                  ParallelDescriptor::MyProc(),
                                  ParallelDescriptor::Communicator());

                if (rc != MPI_SUCCESS)
                    ParallelDescriptor::Abort(rc);
            }
#endif


            Real vol = 1.0;
            for (int dir = 0; dir < BL_SPACEDIM; dir++)
                vol *= amrDataI.DxLevel()[iLevel][dir];

            if (ParallelDescriptor::IOProcessor())
            {
                for (int iComp = 0; iComp < nComp; iComp++)
                {
                    if (norm != 0)
                    {
                        norms[iComp] = norms[iComp] * vol;
                        norms[iComp] = pow(norms[iComp], (1.0/norm));
                    }

                    std::cout << norms[iComp] << " ";
                }
                std::cout << std::endl;
            }
        }

    }
    amrex::Finalize();
}

