clawpack/applications/euler_chem/1d/ModelDetonation/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 "eulerrhok1.h"

#define NEQUATIONS  10
#define NEQUSED     7
#define NFIXUP      6
#define NAUX        0

#include "ClpProblem.h"

#define OWN_FLAGGING
#define OWN_AMRSOLVER
#include "ClpStdProblem.h"
#include "F77Interfaces/F77UpdateLevelTransfer.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 VectorType::InternalDataType DataType;
  typedef AMRSolver<VectorType,FixupType,FlagType,DIM> base;
  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), OutputPMax(0), PThreshold(0.) {
      SetLevelTransfer(new F77UpdateLevelTransfer<VectorType,DIM>(f_prolong, f_restrict, f_tupdate));
    SetFileOutput(new output_type(f_out)); 
    SetFixup(new FixupSpecific(integ));
    SetFlagging(new FlaggingSpecific(*this)); 
  }  

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

  virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
    base::register_at(Ctrl, prefix);
    RegisterAt(LocCtrl,"OutputPMax",OutputPMax); 
    RegisterAt(LocCtrl,"PThreshold",PThreshold); 
  }
  virtual void register_at(ControlDevice& Ctrl) {
    register_at(Ctrl, "");
  }

  virtual void AfterLevelStep(const int Level) {
    if (Level<FineLevel(base::GH())) return;
    if (OutputPMax) {
      double maxp[2];
      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]);
        maxp[0] = 0.0; maxp[1] = 0.0; 
        forall (base::Work(),Time,l,c)       
          BBox bb = base::Work()(Time,l,c).bbox();
          BeginFastIndex1(work, bb, base::Work()(Time,l,c).data(), DataType);
          for (int n=bb.upper(0)-bb.stepsize(0); n>=bb.lower(0)+bb.stepsize(0); n-=bb.stepsize(0)) 
            if (FastIndex1(work,n)>FastIndex1(work,n+bb.stepsize(0)) &&
                FastIndex1(work,n)>FastIndex1(work,n-bb.stepsize(0)) &&
                FastIndex1(work,n)>PThreshold) {
              if (n>maxp[0]) { maxp[0] = double(n); maxp[1] = FastIndex1(work,n); }
              break;
            }
          EndFastIndex1(work);
        end_forall 
      }
#ifdef DAGH_NO_MPI
#else
      if (comm_service::dce() && comm_service::proc_num() > 1) {        
        int num = comm_service::proc_num();       
        int R;
        int size = 2*sizeof(double);
        void *values = (void *) new char[size*num];     
        R = MPI_Allgather(&maxp, size, MPI_BYTE, values, size, MPI_BYTE, 
                          comm_service::comm());
        if ( MPI_SUCCESS != R ) 
          comm_service::error_die( "SolverControlSpecific::maxp", "MPI_Allgather", R );
        for (int i=0;i<num;i++) {
          double tmaxn = *((double *)values + 2*i);
          double tmaxp = *((double *)values + 2*i+1);
          if (tmaxn>maxp[0]) maxp[1] = tmaxp;
        }
        if (values) delete [] (char *)values;
      }
#endif
      int me = MY_PROC; 
      if (me == VizServer) {
        std::ofstream outfile;
        std::ostream* out;
        std::string fname = "pmax.dat";
        outfile.open(fname.c_str(), std::ios::out | std::ios::app);
        out = new std::ostream(outfile.rdbuf());
        *out << base::t[Level] << " " << maxp[1] << std::endl;   
        outfile.close();
        delete out;
      }      
    }    
  }

protected:
  int OutputPMax;
  double PThreshold;
};