clawpack/applications/euler_chem/2d/CellStrucModelDet/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 10
#define NEQUSED     9
#define NFIXUP      8
#define NAUX        4

#define GLOBALMAXIMA  5

#include "ClpProblem.h"

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

#define f_rcveloc FORTRAN_NAME(rcveloc, RCVELOC)
extern "C" {
  void f_rcveloc(const DOUBLE v[]);
}

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 AMRGridAccumulation<DataType,DIM> accumulation_type;
  typedef F77FileOutput<VectorType,DIM> output_type;
public:
  typedef base::vec_grid_data_type vec_grid_data_type;  
  typedef base::grid_data_type grid_data_type;

  SolverSpecific(IntegratorSpecific& integ, 
                 base::initial_condition_type& init,
                 base::boundary_conditions_type& bc) : base(integ, init, bc) {
    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)); 
    _Accumulation = new accumulation_type(f_rcveloc);
    MaxType = 1;
    TrackLevel = -1;
    TrackMin = 0.1;
    TrackMax = 1.e4;
    TrackOutputSingleFile = 1;
    TrackOutputEvery = 1;
    for (register int i=0; i<GLOBALMAXIMA; i++) 
      TrackMaxima[i] = -1;
  }  

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

  virtual void register_at(ControlDevice& Ctrl) { base::register_at(Ctrl); }
  virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
    base::register_at(Ctrl, prefix);
    TrackCtrl = base::LocCtrl.getSubDevice("TrackFront");
    RegisterAt(TrackCtrl,"Level",TrackLevel);
    RegisterAt(TrackCtrl,"Min",TrackMin);
    RegisterAt(TrackCtrl,"Max",TrackMax);
    RegisterAt(TrackCtrl,"OutputSingleFile",TrackOutputSingleFile);
    RegisterAt(TrackCtrl,"OutputEvery",TrackOutputEvery);
    RegisterAt(base::LocCtrl,"MaxType",MaxType);
    char VariableName[32];
    for (register int i=0; i<GLOBALMAXIMA; i++) {
      std::sprintf(VariableName,"Maximum(%d)",i+1);
      RegisterAt(TrackCtrl,VariableName,TrackMaxima[i]);
    }
    _Accumulation->register_at(base::LocCtrl, prefix);
  }

  virtual void SetupData() {
    base::SetupData();
    _Accumulation->SetupData(base::PGH(), base::NGhosts());
    if (TrackLevel>MaxLevel(base::GH())) 
      TrackLevel=MaxLevel(base::GH());
  }

  virtual void Initialize_(const double& dt_start) {
    base::Initialize_(dt_start);
    _Accumulation->Initialize();
  }

  virtual void Output() {
    int Time = CurrentTime(base::GH(),0);
    if (Time == base::LastOutputTime()) return;
    _Accumulation->GridCombine(VizServer);
    int me = MY_PROC; 
    if (me == VizServer) {
      int Level = _Accumulation->AccLevel();
      base::FileOutput_().WriteOut(_Accumulation->AccGrid(),_Accumulation->AccName(),
                                   Time,Level,base::t[0]);
    }
    base::Output();
  }

  virtual void Checkpointing_(const char* CheckpointFile) {
    base::Checkpointing_(CheckpointFile);
    _Accumulation->GridCombine(VizServer);
    if (MY_PROC == VizServer) 
      _Accumulation->Checkpointing();
  }

  virtual bool Restart_(const char* CheckpointFile) {
    if (base::Restart_(CheckpointFile)) {
      if (MY_PROC == VizServer) 
        _Accumulation->Restart();
      else 
        _Accumulation->Initialize();
      return true;
    }
    else 
      return false;
  }

  inline DataType sgn(const DataType& v) const {
    return (v>0. ? 1. : -1.);
  }

  virtual void AdvanceLevel(const int Level, int RegridEvery, bool RegridDone, 
                            bool ShadowAllowed, bool DoFixup,
                            bool RecomposeBaseLev, bool RecomposeHighLev) {

    base::AdvanceLevel(Level,RegridEvery,RegridDone,ShadowAllowed,
                       DoFixup,RecomposeBaseLev,RecomposeHighLev);

    int Time = CurrentTime(base::GH(),Level);
    int TStep = TimeStep(base::U(),Level);
    if (MaxType) {
      AdaptiveBoundaryUpdate(base::U(),Time,Level);
      Sync(base::U(),Time,Level);
      ExternalBoundaryUpdate(base::U(),Time,Level);
      
      int irho = 1;
      int iu   = NEQUSED-DIM-2;
      int iv   = NEQUSED-DIM-1;
      ((output_type*) base::_FileOutput)->Transform(base::U(), base::Work(), Time+TStep, 
                                                    Level, irho, base::t[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);
        BeginFastIndex2(WorkN, Work()(Time+TStep,Level,c).bbox(), 
                        Work()(Time+TStep,Level,c).data(), DataType);
        DataType uy, vx;
        DCoords dx = base::GH().worldStep(ss);
        for_2 (i, j, bbi, ss)
          vx = (FastIndex2(U,i+ss(0),j)(iv)/FastIndex2(WorkN,i+ss(0),j)-
                FastIndex2(U,i-ss(0),j)(iv)/FastIndex2(WorkN,i-ss(0),j))/(2.0*dx(0));
          uy = (FastIndex2(U,i,j+ss(1))(iu)/FastIndex2(WorkN,i,j+ss(1))-
                FastIndex2(U,i,j-ss(1))(iu)/FastIndex2(WorkN,i,j-ss(1)))/(2.0*dx(1));
          FastIndex2(Work,i,j) = std::fabs(vx-uy);
        end_for
        EndFastIndex2(U);      
        EndFastIndex2(Work);      
        EndFastIndex2(WorkN);      
      end_forall
    }
    else {
      int ipress = base::Dim()+4;      
      ((output_type*) base::_FileOutput)->Transform(base::U(), base::Work(), Time, Level, 
                                                    ipress, base::t[Level]);
    }
  
    _Accumulation->Accumulate(base::Work(), Time,  Level, base::t[Level]);
  }

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

    double cfl = base::Tick(VariableTimeStepping,dtv,cflv,Rejections);
    if (TrackLevel<0) return cfl;

    int TimeZero  = CurrentTime(base::GH(),0);
    TimeZero /= RefinedBy(base::GH(),MaxLevel(base::GH()));
    if (TimeZero % TrackOutputEvery != 0) return cfl;

    int Level = std::min(TrackLevel,FineLevel(base::GH()));
    int Time = CurrentTime(base::GH(),Level);
    int TStep = TimeStep(base::U(),Level);

    BBox wholebb(base::GH().wholebbox());
    wholebb.refine(base::GH().refinedby(Level));
    int idxl = wholebb.extents(1);
    
    {
      // Track the shock front in each row of cells as greatest x-location 
      // with local minimum in curvature. This captures basically the von
      // Neumann point.
      int fds = base::Dim();
      DataType* front_data = new DataType[(base::NEquations()+fds)*idxl];
      VectorType* front_val = (VectorType *)(front_data+fds*idxl);
      for (register int j=0;j<idxl; j++) {
        for (register int i=0;i<fds; i++) 
          front_data[fds*j+i]   = -1.e37;
        front_val[j] = -1.e37;
      }
      
      int ipress = base::Dim()+4;      
      ((output_type*) base::_FileOutput)->Transform(base::U(), base::Work(), Time, Level, 
                                                    ipress, base::t[Level]);
      forall (base::Work(),Time,Level,c)
        Coords ss(base::Work()(Time,Level,c).bbox().stepsize());
        BBox bbi(base::Work()(Time,Level,c).bbox());
        bbi.upper(0) -= bbi.stepsize(0);
        BeginFastIndex2(P, Work()(Time,Level,c).bbox(), 
                        Work()(Time,Level,c).data(), DataType);
        BeginFastIndex2(dP, Work()(Time+TStep,Level,c).bbox(), 
                        Work()(Time+TStep,Level,c).data(), DataType);
        DCoords dx = base::GH().worldStep(ss);
        // 1st derivative
        for_2 (i, j, bbi, ss)
          FastIndex2(dP,i,j) = (FastIndex2(P,i+ss(0),j)-FastIndex2(P,i,j))/dx(0);
        end_for
        bbi.lower(0) += bbi.stepsize(0);
        grid_data_type DataddP(bbi);
        BeginFastIndex2(ddP, DataddP.bbox(), DataddP.data(), DataType);
        // 2nd derivative
        for_2 (i, j, bbi, ss)
          FastIndex2(ddP,i,j) = (FastIndex2(dP,i,j)-FastIndex2(dP,i-ss(0),j))/dx(0);
        end_for
        bbi.upper(0) -= bbi.stepsize(0);
        // 3rd derivative
        for_2 (i, j, bbi, ss)
          FastIndex2(dP,i,j) = (FastIndex2(ddP,i+ss(0),j)-FastIndex2(ddP,i,j))/dx(0);
        end_for
        bbi = base::U().interiorbbox(Time,Level,c);
        for_2 (i, j, bbi, ss)
          if (FastIndex2(P,i,j)>=TrackMin && FastIndex2(P,i,j)<=TrackMax &&
              FastIndex2(ddP,i,j)<0. && 
              sgn(FastIndex2(dP,i,j))!=sgn(FastIndex2(dP,i-ss(0),j))) {
            DCoords fc = base::GH().worldCoords(Coords(2,i,j),ss);
            fc(0) *= (1.+dx(0)/(FastIndex2(dP,i-ss(0),j)-FastIndex2(dP,i,j)));
            int jj = (j-wholebb.lower(1))/ss(1);
            if (fc(0) > front_data[2*jj]) {
              front_data[fds*jj] = fc(0);
              front_data[fds*jj+1] = fc(1)+0.5*dx(1);
              front_val[jj] = base::U()(Time,Level,c)(i,j);
              DataType pmax = FastIndex2(P,i,j);
              for (int ii=i-base::NGhosts()*ss(0); ii<=i+base::NGhosts()*ss(0); ii+=ss(0)) {
                if (FastIndex2(P,ii,j) > pmax) {
                  front_val[jj] = base::U()(Time,Level,c)(ii,j);
                  pmax = FastIndex2(P,ii,j);
                }
              }
            }
          }
        end_for
        EndFastIndex2(P);      
        EndFastIndex2(dP);      
        EndFastIndex2(ddP);      
      end_forall

#ifdef DAGH_NO_MPI
#else
      if (comm_service::dce() && comm_service::proc_num() > 1) {        
        int num = comm_service::proc_num();
        MPI_Status status;
        int R;
        int size = sizeof(DataType)*(NEquations()+fds)*idxl;
        int tag = 10000;
        if (MY_PROC!=VizServer) {
          R = MPI_Send((void *)front_data, size, MPI_BYTE, VizServer, tag, comm_service::comm());
          if ( MPI_SUCCESS != R ) 
            comm_service::error_die( "Tick", "MPI_Send", R );
        }
        else 
          for (int proc=0; proc<num; proc++) {
            if (proc==VizServer) continue;
            DataType* dat = new DataType[(NEquations()+fds)*idxl];
            VectorType* datv = (VectorType *)(dat+fds*idxl);
            R = MPI_Recv((void *)dat, size, MPI_BYTE, proc, tag, comm_service::comm(), &status);
            if ( MPI_SUCCESS != R ) 
              comm_service::error_die( "Tick", "MPI_Recv", R );
            for (register int j=0; j<idxl; j++) 
              if (dat[fds*j] > front_data[fds*j]) {
                for (register int i=0;i<fds; i++) 
                  front_data[fds*j+i] = dat[fds*j+i];
                front_val[j] = datv[j];
              }
            delete [] dat;
          }
      }
#endif

      int me = MY_PROC; 
      if (me == VizServer) {
        DataType* output_data = new DataType[_FileOutput->Ncnt()*idxl];
        BBox bb(2,0,0,idxl-1,0,1,1);
        vec_grid_data_type val(bb,front_val);
        for (int cnt=1;cnt<=_FileOutput->Ncnt(); cnt++) {
          if (_FileOutput->OutputName(cnt-1).c_str()[0] == '-') continue;
          grid_data_type output(bb,output_data+(cnt-1)*idxl);
          ((output_type::out_2_func_type) (((output_type*) base::_FileOutput)->GetOutFunc()))
            (FA(2,val),FA(2,output),BOUNDING_BOX(bb),base::NEquations(),cnt,base::t[Level]);
          DataType* dummy;
          output.deallocate(dummy);
        }
        std::ofstream outfile;
        std::ostream* out;
        char name[256];
        if (TrackOutputSingleFile) 
          sprintf(name,"front.txt");
        else 
          sprintf(name,"front_%d.txt",Time);
        std::string fname = name;
        outfile.open(fname.c_str(), std::ios::out | std::ios::app);
        out = new std::ostream(outfile.rdbuf());
        if (TrackOutputSingleFile)        
          *out << "**********   t=" << base::t[Level] << "   Step:" 
               << Time << "   **********" << std::endl;
        for (register int j=0;j<idxl; j++) {
          *out << front_data[fds*j] << " " << front_data[fds*j+1];
          for (int cnt=0;cnt<_FileOutput->Ncnt(); cnt++) 
            if (_FileOutput->OutputName(cnt).c_str()[0] != '-') 
              *out << " " << output_data[cnt*idxl+j];
          *out << std::endl;
        }
        val.deallocate(front_val);
        delete [] output_data;
        outfile.close();
        delete out;      
      }
    }

    // Track the maximum of selectable output components in the entire domain.
    for (register int im=0; im<GLOBALMAXIMA; im++) 
      if (TrackMaxima[im]>=0) {
        int fds = base::Dim()+1;
        DataType* front_data = new DataType[(base::NEquations()+fds)*idxl];
        VectorType* front_val = (VectorType *)(front_data+fds*idxl);
        for (register int j=0;j<idxl; j++) {
          for (register int i=0;i<fds; i++) 
          front_data[fds*j+i]   = -1.e37;
          front_val[j] = -1.e37;
        }
        
        ((output_type*) base::_FileOutput)->Transform(base::U(), base::Work(), Time, Level, 
                                                      TrackMaxima[im], base::t[Level]);
        forall (base::Work(),Time,Level,c)
          Coords ss(base::Work()(Time,Level,c).bbox().stepsize());
          BBox bbi(base::Work()(Time,Level,c).bbox());
          bbi.upper(0) -= bbi.stepsize(0);
          BeginFastIndex2(V, Work()(Time,Level,c).bbox(), 
                          Work()(Time,Level,c).data(), DataType);
          DCoords dx = base::GH().worldStep(ss);
          bbi = base::U().interiorbbox(Time,Level,c);
          for_2 (i, j, bbi, ss)
            int jj = (j-wholebb.lower(1))/ss(1);
            if (FastIndex2(V,i,j)>front_data[fds*jj]) {
              front_data[fds*jj] = FastIndex2(V,i,j);
              DCoords fc = base::GH().worldCoords(Coords(2,i,j),ss);
              front_data[fds*jj+1] = fc(0)+0.5*dx(0);
              front_data[fds*jj+2] = fc(1)+0.5*dx(1);
              front_val[jj] = base::U()(Time,Level,c)(i,j);
            }
          end_for
          EndFastIndex2(V);      
        end_forall

#ifdef DAGH_NO_MPI
#else
        if (comm_service::dce() && comm_service::proc_num() > 1) {      
          int num = comm_service::proc_num();
          MPI_Status status;
          int R;
          int size = sizeof(DataType)*(NEquations()+fds)*idxl;
          int tag = 10000;
          if (MY_PROC!=VizServer) {
            R = MPI_Send((void *)front_data, size, MPI_BYTE, VizServer, tag, comm_service::comm());
            if ( MPI_SUCCESS != R ) 
              comm_service::error_die( "Tick", "MPI_Send", R );
          }
          else 
            for (int proc=0; proc<num; proc++) {
              if (proc==VizServer) continue;
              DataType* dat = new DataType[(NEquations()+fds)*idxl];
              VectorType* datv = (VectorType *)(dat+fds*idxl);
              R = MPI_Recv((void *)dat, size, MPI_BYTE, proc, tag, comm_service::comm(), &status);
              if ( MPI_SUCCESS != R ) 
                comm_service::error_die( "Tick", "MPI_Recv", R );
              for (register int j=0; j<idxl; j++) 
                if (dat[fds*j] > front_data[fds*j]) {
                  for (register int i=0;i<fds; i++) 
                    front_data[fds*j+i] = dat[fds*j+i];
                  front_val[j] = datv[j];
                }
              delete [] dat;
            }
        }
#endif

        int me = MY_PROC; 
        if (me == VizServer) {
          DataType* output_data = new DataType[_FileOutput->Ncnt()*idxl];
          BBox bb(2,0,0,idxl-1,0,1,1);
          vec_grid_data_type val(bb,front_val);
          for (int cnt=1;cnt<=_FileOutput->Ncnt(); cnt++) {
            if (_FileOutput->OutputName(cnt-1).c_str()[0] == '-') continue;
            grid_data_type output(bb,output_data+(cnt-1)*idxl);
            ((output_type::out_2_func_type) (((output_type*) base::_FileOutput)->GetOutFunc()))
              (FA(2,val),FA(2,output),BOUNDING_BOX(bb),base::NEquations(),cnt,base::t[Level]);
            DataType* dummy;
            output.deallocate(dummy);
          }
          std::ofstream outfile;
          std::ostream* out;
          char name[256];
          if (TrackOutputSingleFile) 
            sprintf(name,"max_%d.txt",TrackMaxima[im]);
          else 
            sprintf(name,"max_%d_%d.txt",TrackMaxima[im],Time);
          std::string fname = name;
          outfile.open(fname.c_str(), std::ios::out | std::ios::app);
          out = new std::ostream(outfile.rdbuf());
          if (TrackOutputSingleFile)        
            *out << "**********   t=" << base::t[Level] << "   Step:" 
                 << Time << "   **********" << std::endl;
          for (register int j=0;j<idxl; j++) {
          *out << front_data[fds*j+1] << " " << front_data[fds*j+2];
          for (int cnt=0;cnt<_FileOutput->Ncnt(); cnt++) 
            if (_FileOutput->OutputName(cnt).c_str()[0] != '-') 
              *out << " " << output_data[cnt*idxl+j];
          *out << std::endl;
          }
          val.deallocate(front_val);
          delete [] output_data;
          outfile.close();
          delete out;      
        }
        delete [] front_data;
      }

    return cfl;
  }

protected:
  int MaxType, TrackLevel, TrackOutputEvery, TrackOutputSingleFile;
  int TrackMaxima[GLOBALMAXIMA];
  double TrackMin, TrackMax;
  accumulation_type* _Accumulation;
  ControlDevice TrackCtrl;