clawpack/applications/intro/3d/Refinement/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 "DAGH.h"                           
#include "IO/control-device.h"
#include "Vector.h"
#include "VectorLarge.h"

#define DIM        3
#define NEQUATIONS 1
typedef DOUBLE DataType;
typedef Vector<DataType,NEQUATIONS> VectorType;

#include "DAGHFortranInterfaces.h"                           

//----------------------------------------------------------------------
// Dummy integrator
//----------------------------------------------------------------------

#include "Integrator.h"  

class IntegratorSpecific : public Integrator<VectorType,DIM> {
  typedef Integrator<VectorType,DIM> base;
public:
  IntegratorSpecific() : _MethodOrder(2) {}

  virtual double CalculateGrid(base::vec_grid_data_type& NewStateVec, 
                               base::vec_grid_data_type& OldStateVec,
                               base::vec_grid_data_type* Flux[],
                               const int& level, 
                               const double& t, const double& dt, 
                               const int& mpass) 
  { NewStateVec.copy(OldStateVec); return 1.0; } 

  virtual void AllocGridFluxes(const BBox &bb, base::vec_grid_data_type**& Flux) {
    Flux = new base::vec_grid_data_type* [2*base::Dim()]; 
    for (register int d=0; d<base::Dim(); d++) { 
      Flux[2*d] = new base::vec_grid_data_type(bb); 
      Flux[2*d+1] = Flux[2*d];
    }
  }

  virtual void DeAllocGridFluxes(base::vec_grid_data_type**& Flux) {
    if (Flux) {
      for (register int d=0; d<base::Dim(); d++) 
        if (Flux[2*d]) delete Flux[2*d];
      delete [] Flux;
    }
  }
    
  virtual void ResetGridFluxes(vec_grid_data_type**& Flux) {
    if (Flux) {
      VectorType zero(0.0);
      for (register int d=0; d<2*base::Dim(); d++) 
        if (Flux[d]) 
          Flux[d]->equals(zero);
    }
  }

  virtual int ControlGrid(base::vec_grid_data_type& StateVec, const int& level, 
                          const BBox& where, const double& time, const int verbose) 
        { return 1; }  
  virtual int NMethodOrder() const { return _MethodOrder; }
private:
  int _MethodOrder;
};


//----------------------------------------------------------------------
// Initial conditions in Fortran 77
//----------------------------------------------------------------------

#define f_initial FORTRAN_NAME(ic, IC)

extern "C" {
  void f_initial();
}

#include "F77Interfaces/F77InitialCondition.h"  

class InitialConditionSpecific : 
  public F77InitialCondition<VectorType,DIM> {
public:    
  InitialConditionSpecific() : 
    F77InitialCondition<VectorType,DIM>(f_initial) {}
};

//----------------------------------------------------------------------
// Boundary conditions in Fortran 77
//----------------------------------------------------------------------

#define f_boundary FORTRAN_NAME(physbd, PHYSBD)

extern "C" {
  void f_boundary();
}

#include "F77Interfaces/F77BoundaryConditions.h"  

class BoundaryConditionsSpecific : 
  public F77BoundaryConditions<VectorType,DIM> {
public:    
  BoundaryConditionsSpecific() : 
    F77BoundaryConditions<VectorType,DIM>(f_boundary) {}
};

//----------------------------------------------------------------------
// Inter-level operators in Fortran 77
//----------------------------------------------------------------------

#define f_restrict FORTRAN_NAME(restrict3, RESTRICT3)
#define f_prolong FORTRAN_NAME(prolong3, PROLONG3)

extern "C" {
  void f_restrict();
  void f_prolong();
}

#include "F77Interfaces/F77LevelTransfer.h"  

class LevelTransferSpecific : 
  public F77LevelTransfer<VectorType,DIM> {
public:    
  LevelTransferSpecific() : 
    F77LevelTransfer<VectorType,DIM>(f_prolong, f_restrict) {}
};

//----------------------------------------------------------------------
// Flux correction
//----------------------------------------------------------------------

#define NFIXUP     5   
typedef Vector<DataType,NFIXUP> FixupType;

#include "AMRFixup.h"  

class FixupSpecific : public AMRFixup<VectorType,FixupType,DIM> {
public:
  FixupSpecific() : AMRFixup<VectorType,FixupType,DIM>() {}
};

//----------------------------------------------------------------------
// Flagging for refinement
//----------------------------------------------------------------------

typedef short FlagType; 

#include "AMRFlagging.h"
#include "Criteria/LimiterType.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 LimiterType<VectorType,FlagType,DIM>());
  }

  ~FlaggingSpecific() { DeleteAllCriterions(); }
};

//----------------------------------------------------------------------
// Berger-Collela AMR algorithmus
//----------------------------------------------------------------------

#include "AMRSolver.h" 

class SolverSpecific : 
  public AMRSolver<VectorType,FixupType,FlagType,DIM> {
  typedef AMRSolver<VectorType,FixupType,FlagType,DIM> base;
public:
  SolverSpecific(IntegratorSpecific& integ, 
                 base::initial_condition_type& init,
                 base::boundary_conditions_type& bc) :
    AMRSolver<VectorType,FixupType,FlagType,DIM>(integ, init, bc) {
    SetLevelTransfer(new LevelTransferSpecific());
    SetFileOutput(new FileOutput<VectorType,DIM>()); 
    SetFixup(new FixupSpecific());
    SetFlagging(new FlaggingSpecific(*this)); 
  }  

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

  virtual void Initialize(double& t, double& dt) {  
    // Must be turned on with periodic boundary conditions!
    U().GF_SetAdaptBoundaryType(DAGHNoAdaptBoundary);
    U().GF_SetCommType(DAGHNoComm);
    // U().GF_SetInitializeFlag(DAGHFalse);
    AMRSolver<VectorType,FixupType,FlagType,DIM>::Initialize(t,dt);
  }
}; 
 
//----------------------------------------------------------------------
// Time stepping and file output
//----------------------------------------------------------------------

#include "SolverControl.h"

class SolverControlSpecific : 
  public SolverControl {
public:
  SolverControlSpecific(AMRSolver<VectorType,FixupType,FlagType,DIM>& solver) : 
    SolverControl(solver) {}
};