clawpack/applications/euler_chem/2d/PipeBlast/src/Problem.h

// -*- C++ -*-

// Copyright (C) 2003-2007 California Institute of Technology
// Ralf Deiterding, ralf@amroc.net

#ifndef AMROC_PROBLEM_H
#define AMROC_PROBLEM_H

#include "eulerrhok2.h"
#define NEQUATIONS  6
#define NEQUSED     6
#define NFIXUP      5
#define NAUX        1

#include "ClpProblem.h"

#define OWN_FLAGGING
#define OWN_AMRSOLVER
#include "ClpStdProblem.h"
#include "F77Interfaces/F77UpdateLevelTransfer.h"
#include "AMRInterpolation.h"

class FlaggingSpecific : 
  public AMRFlagging<VectorType,FixupType,FlagType,DIM> {
  typedef AMRFlagging<VectorType,FixupType,FlagType,DIM> base;
public:
  FlaggingSpecific(base::solver_type& solver) : base(solver) {
      base::AddCriterion(new ScaledGradient<VectorType,FlagType,DIM>());
      base::AddCriterion(new LimiterType<VectorType,FlagType,DIM>());
      base::AddCriterion(new AbsoluteError<VectorType,FixupType,FlagType,DIM>(solver));
      base::AddCriterion(new RelativeError<VectorType,FixupType,FlagType,DIM>(solver));
      base::AddCriterion(new F77ScaledGradient<VectorType,FlagType,DIM>(f_flg));
      base::AddCriterion(new F77LimiterType<VectorType,FlagType,DIM>(f_flg));
      base::AddCriterion(new F77AbsoluteError<VectorType,FixupType,FlagType,DIM>(solver,f_flg));
      base::AddCriterion(new F77RelativeError<VectorType,FixupType,FlagType,DIM>(solver,f_flg));
    }
  ~FlaggingSpecific() { DeleteAllCriterions(); }
};

class SolverSpecific : 
  public AMRSolver<VectorType,FixupType,FlagType,DIM> {
  typedef AMRSolver<VectorType,FixupType,FlagType,DIM> base;
  typedef VectorType::InternalDataType DataType;
  typedef AMRInterpolation<VectorType,DIM> interpolation_type;
  typedef interpolation_type::point_type point_type;
  typedef F77FileOutput<VectorType,DIM> output_type;
public:
  SolverSpecific(IntegratorSpecific& integ, 
                 base::initial_condition_type& init,
                 base::boundary_conditions_type& bc) :
    AMRSolver<VectorType,FixupType,FlagType,DIM>(integ, init, bc), OutputPressureEvery(1) {
    SetLevelTransfer(new F77UpdateLevelTransfer<VectorType,DIM>(f_prolong, f_restrict, f_tupdate));
    SetFileOutput(new F77FileOutput<VectorType,DIM>(f_out)); 
    SetFixup(new FixupSpecific(integ));
    SetFlagging(new FlaggingSpecific(*this)); 
    _Interpolation = new interpolation_type();
  }  

  ~SolverSpecific() {
    delete _LevelTransfer;
    delete _Flagging;
    delete _Fixup;
    delete _FileOutput;
    delete _Interpolation;
  }

  virtual void SetupData() {
    base::SetupData();
    _Interpolation->SetupData(base::PGH(), base::NGhosts());
  }
  virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
    base::register_at(Ctrl,prefix);
    RegisterAt(base::LocCtrl,"OutputPressureEvery",OutputPressureEvery);
  } 
  virtual void register_at(ControlDevice& Ctrl) {
    base::register_at(Ctrl);
  }
  // Overloading of Tick() to evaluate force integral on spheres after each
  // level 0 time step.

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

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

    int TimeZero  = CurrentTime(base::GH(),0);
    TimeZero /= RefinedBy(base::GH(),MaxLevel(base::GH()));
    if (TimeZero%OutputPressureEvery!=0 || 
        base::GH().nbndry()<=1) return cfl;

    int me = MY_PROC; 

    // Calculate pressure
    for (register int l=0; l<=FineLevel(base::GH()); l++) {
      int Time = CurrentTime(base::GH(),l);
      int press_cnt = base::Dim()+4;      
      ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time, l, 
                                              press_cnt, base::t[l]);
    }

    const double* bndrycoords = base::GH().wholebndry();
    double lpy = bndrycoords[1+base::GH().nbndry()*(0+2*1)];
    double upy = bndrycoords[1+base::GH().nbndry()*(1+2*1)];
    double upx = bndrycoords[1+base::GH().nbndry()*(1+2*0)];
    double dy = base::GH().delta_x(1,MaxLevel(base::GH()));
    double pav = 0., global_pav = 0.;
    int Npoints = int((upy-lpy-0.5*dy)/dy);

    DataType* p = new DataType[Npoints];
    point_type* xc = new point_type[Npoints];
    register int n;
    for (n=0; n<Npoints; n++) {
      xc[n](0) = upx;
      xc[n](1) = lpy+(n+0.5)*dy;
    }
    
    _Interpolation->PointsValues(base::Work(),base::t[0],Npoints,xc,p);
    for (n=0; n<Npoints; n++)
      if (p[n]>-1.e37) 
        pav += p[n]; 
    delete [] p;
    delete [] xc;

#ifdef DAGH_NO_MPI
    global_pav = pav;
#else
    int Nval = 1;
    MPI_Allreduce(&pav, &global_pav, Nval, MPI_DOUBLE, MPI_SUM, base::Comm());
#endif

    // Append to output file only on processor VizServer
    if (me == VizServer) {
      std::ofstream outfile;
      std::ostream* out;
      std::string fname = "pressure.dat";
      outfile.open(fname.c_str(), std::ios::out | std::ios::app);
      out = new std::ostream(outfile.rdbuf());
      *out << base::t[0] << " " << global_pav/Npoints << std::endl;
      outfile.close();
      delete out;      
    }

    return cfl;
  } 

protected:
  int OutputPressureEvery;
  interpolation_type* _Interpolation; 
};