weno/applications/euler/2d/Vortex/src/Problem.h

// -*- C++ -*-

#ifndef AMROC_PROBLEM_H
#define AMROC_PROBLEM_H

#define DIM  2

//#define SYNCTIME    
#define OWN_AMRSOLVER
#include "WENOProblem.h"

#define f_exact FORTRAN_NAME(exact, EXACT)

extern "C" {
  void f_exact(); 
}

#ifdef SYNCTIME
#include "WENOStdSyncTimeProblem.h"
#else
#include "WENOStdProblem.h"
#endif
#include "F77Interfaces/F77ExactSolution.h"

class SolverSpecific : 
#ifdef SYNCTIME
  public AMRPreAdaptSolverSyncTime<VectorType,FixupType,FlagType,DIM> {
  typedef AMRPreAdaptSolverSyncTime<VectorType,FixupType,FlagType,DIM> base;
#else
  public AMRPreAdaptSolver<VectorType,FixupType,FlagType,DIM> {
  typedef AMRPreAdaptSolver<VectorType,FixupType,FlagType,DIM> base;
#endif
  typedef AMRSolver<VectorType,FixupType,FlagType,DIM> bbase;
  typedef WENOIntegrator<VectorType,DIM> weno_integ_type;
  typedef VectorType::InternalDataType DataType;
public:
  SolverSpecific(IntegratorSpecific& integ, 
                 bbase::initial_condition_type& init,
                 bbase::boundary_conditions_type& bc) :
#ifdef SYNCTIME
    AMRPreAdaptSolverSyncTime<VectorType,FixupType,FlagType,DIM>(integ, init, bc) {
#else
    AMRPreAdaptSolver<VectorType,FixupType,FlagType,DIM>(integ, init, bc) {
#endif
    SetLevelTransfer(new F77LevelTransfer<VectorType,DIM>(f_prolong, f_restrict));
#ifdef f_flgout
    SetFileOutput(new F77FileOutput<VectorType,DIM>(f_flgout)); 
#else   
    SetFileOutput(new FileOutput<VectorType,DIM>()); 
#endif
    SetFixup(new FixupSpecific());
    SetFlagging(new FlaggingSpecific(*this)); 

    SetExactSolution(new F77ExactSolution<VectorType,DIM>(f_exact)); 
  }

  ~SolverSpecific() {
    delete _LevelTransfer;
    delete _Flagging;
    delete _Fixup;
    delete _FileOutput;
  }
  
  virtual void Initialize_(const double& dt_start) {   
    base::Initialize_(dt_start);
    CalculateCurrentMass(mass_old);
  }

  virtual double Tick(int VariableTimeStepping, const double dtv[], 
                      const double cflv[], int& Rejections) {

    double cfl = base::Tick(VariableTimeStepping,dtv,cflv,Rejections);

    VectorType mass_new;
    CalculateCurrentMass(mass_new);

    // Append to output file only on processor VizServer
    int me = MY_PROC; 
    if (me == VizServer) {
      std::ofstream outfile;
      std::ostream* out;
      std::string fname = "mass.dat";
      outfile.open(fname.c_str(), std::ios::out | std::ios::app);
      out = new std::ostream(outfile.rdbuf());
      *out << base::t[0];
      for (int n=0; n<base::NEquations(); n++)
        *out << " " << mass_old[n]-mass_new[n];
      *out << std::endl;
      outfile.close();
      delete out;      
    }

    // compute kinetic energy
    int irho = 0;
    int iu   = 1;
    int iv   = 2;
    double kin = 0.0;
    for ( int Level = FineLevel(base::GH()) ; Level >= 0 ; Level --)
      {
        int Time = CurrentTime(base::GH(),Level);
        forall (base::U(),Time,Level,c)
           Coords ss(base::U()(Time,Level,c).bbox().stepsize());
           BBox bbi(base::U().interiorbbox(Time,Level,c));
           BeginFastIndex2(U, base::U()(Time,Level,c).bbox(), 
                           base::U()(Time,Level,c).data(), VectorType);
           BeginFastIndex2(Work, Work()(Time,Level,c).bbox(), 
                           Work()(Time,Level,c).data(), DataType);
           DataType rho, u, v;
           for_2 (i, j, bbi, ss)
              rho = FastIndex2(U,i,j)(irho);
              u = FastIndex2(U,i,j)(iu)/rho;
              v = FastIndex2(U,i,j)(iv)/rho;
              FastIndex2(Work,i,j) = (u*u+v*v)/2.0;
           end_for
           EndFastIndex2(U);
           EndFastIndex2(Work);
        end_forall

        forall (base::Work(),Time,Level,c) 
             if (Level<FineLevel(base::GH())) {
                forall (base::Work(),Time,Level+1,c2) 
                    BBox bb(coarsen(base::Work().interiorbbox(Time,Level+1,c2),base::GH().refinefactor(Level)));
                    base::Work()(Time,Level,c).equals(DataType(0.0),bb);
                end_forall      
             }
        end_forall      
        kin += base::Work().GF_sum(Time,Level)/(Level+1)/(Level+1);
      }
    DCoords dx = base::GH().worldStep(base::U().GF_StepSize(0));
    for (int d=0; d<base::Dim(); d++) kin *= dx(d);

    if (MY_PROC == VizServer) {
      std::ofstream outfile;
      std::ostream* out;
      std::string fname = "stats.dat";
      outfile.open(fname.c_str(), std::ios::out | std::ios::app);
      out = new std::ostream(outfile.rdbuf());
      out->setf(std::ios::scientific, std::ios::floatfield);
      *out << base::t[0] << " " << kin << std::endl;
      outfile.close();
      delete out;      
    }

    return cfl;
  }

  void CalculateCurrentMass(VectorType& mass) {
    int Level=0;
    int Time = CurrentTime(base::GH(),Level); 
    int TStep = TimeStep(base::U(),Level);      
    forall (base::U(),Time,Level,c)   
      base::U()(Time+TStep,Level,c).equals(base::U()(Time,Level,c));
    end_forall

    DCoords dx = base::GH().worldStep(base::U().GF_StepSize(Level));
    mass = Sum(base::U(),Time+TStep,Level);
    for (int d=0; d<base::Dim(); d++) mass *= dx(d);
  }

protected:
VectorType mass_old; 
};