weno/applications/euler/2d/TestCBC/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;      
    }

    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; 
};