clawpack/applications/euler_znd/2d/StatDet/src/Problem.h

// -*- C++ -*-

// Copyright (C) 2002 Ralf Deiterding
// Brandenburgische Universitaet Cottbus

#ifndef AMROC_PROBLEM_H
#define AMROC_PROBLEM_H

#include "eulerznd2.h"

#undef  NAUX
#define NAUX      4

#include "ClpProblem.h"

#define OWN_AMRSOLVER
#include "ClpStdProblem.h"
#include "SpecialOutput/AMRGridAccumulation.h"

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

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:
  SolverSpecific(IntegratorSpecific& integ, 
                 base::initial_condition_type& init,
                 base::boundary_conditions_type& bc) :
    AMRSolver<VectorType,FixupType,FlagType,DIM>(integ, init, bc) {
    SetLevelTransfer(new F77LevelTransfer<VectorType,DIM>(f_prolong, f_restrict));
    SetFileOutput(new F77FileOutput<VectorType,DIM>(f_flgout)); 
    SetFixup(new FixupSpecific(integ));
    SetFlagging(new FlaggingSpecific(*this)); 
    _Accumulation = new accumulation_type(f_rcveloc);
    MaxType = 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);
    RegisterAt(base::LocCtrl,"MaxType",MaxType);
    _Accumulation->register_at(base::LocCtrl, prefix);
  }

  virtual void SetupData() {
    base::SetupData();
    _Accumulation->SetupData(base::PGH(), base::NGhosts());
  }

  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->Checkpointing();
  }

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

  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);
    if (MaxType) {
      AdaptiveBoundaryUpdate(base::U(),Time,Level);
      Sync(base::U(),Time,Level);
      ExternalBoundaryUpdate(base::U(),Time,Level);
      
      int irho = 1;
      int iu   = 2;
      int iv   = 3;
      int TStep = TimeStep(base::U(),Level);
      forall (base::U(),Time,Level,c)
        ((output_type*) base::_FileOutput)->Transform(base::U(), base::Work(), Time+TStep, 
                                                      Level, irho, base::t[Level]);
        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 {
      forall (base::U(),Time,Level,c)
        int ipress = base::Dim()+4;      
        ((output_type*) base::_FileOutput)->Transform(base::U(), base::Work(), Time, Level, 
                                                      ipress, base::t[Level]);
      end_forall
    }
  
    _Accumulation->Accumulate(base::Work(), Time,  Level, base::t[Level]);
  }
private:
  int MaxType;
  accumulation_type* _Accumulation;