clawpack/applications/euler/2d/Spheres/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 <cmath>
#include "euler2.h"
#include "ClpProblem.h"
#undef NAUX
#define NAUX 4

#define f_combl_read FORTRAN_NAME(cblread, CBLREAD)

extern "C" {
  void f_combl_read(INTEGER& NSph, DOUBLE xm[], DOUBLE rd[], INTEGER Np[]); 
}

#define OWN_FLAGGING
#define OWN_GFMAMRSOLVER
#include "ClpStdGFMProblem.h"
#include "AMRGFMInterpolation.h"

class FlaggingSpecific : 
  public AMRFlagging<VectorType,FixupType,FlagType,DIM> {
  typedef AMRFlagging<VectorType,FixupType,FlagType,DIM> base;
public:
  FlaggingSpecific(base::solver_type& solver) : base(solver), _MaxLev(10000), _LastTime(10000.0) {
      base::AddCriterion(new ByValue<VectorType,FlagType,DIM>());
      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));
#ifdef f_flgout
      base::AddCriterion(new F77ByValue<VectorType,FlagType,DIM>(f_flgout));
      base::AddCriterion(new F77ScaledGradient<VectorType,FlagType,DIM>(f_flgout));
      base::AddCriterion(new F77LimiterType<VectorType,FlagType,DIM>(f_flgout));
      base::AddCriterion(new F77AbsoluteError<VectorType,FixupType,FlagType,DIM>(solver,f_flgout));
      base::AddCriterion(new F77RelativeError<VectorType,FixupType,FlagType,DIM>(solver,f_flgout));
#endif
    }

  ~FlaggingSpecific() { DeleteAllCriterions(); }

  virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
    base::register_at(Ctrl, prefix);
    RegisterAt(base::LocCtrl,"MaxLev",_MaxLev);    
    RegisterAt(base::LocCtrl,"LastTime",_LastTime);    
  }
  virtual void register_at(ControlDevice& Ctrl) { register_at(Ctrl, ""); }

  virtual void SetFlags(const int Time, const int Level, double t, double dt) {
    if (t>_LastTime || Level<=_MaxLev) 
       base::SetFlags(Time, Level, t,dt);
    else {
      forall (base::Flags(),Time,Level,c)  
        base::Flags()(Time,Level,c) = GoodPoint;
      end_forall
    }
  }

protected:
  int _MaxLev;
  double _LastTime;
};


class SolverSpecific : 
  public AMRGFMSolver<VectorType,FixupType,FlagType,DIM> {
  typedef VectorType::InternalDataType DataType;
  typedef AMRGFMSolver<VectorType,FixupType,FlagType,DIM> base;
  typedef AMRGFMInterpolation<VectorType,FixupType,FlagType,DIM> interpolation_type;
  typedef interpolation_type::point_type point_type;
  typedef F77GFMFileOutput<VectorType,FixupType,FlagType,DIM> output_type;
public:
  SolverSpecific(IntegratorSpecific& integ, 
                 base::initial_condition_type& init,
                 base::boundary_conditions_type& bc) : base(integ, init, bc), OutputDragEvery(1) {
    SetLevelTransfer(new F77LevelTransfer<VectorType,DIM>(f_prolong, f_restrict));
#ifdef f_flgout
    SetFileOutput(new F77GFMFileOutput<VectorType,FixupType,FlagType,DIM>(*this,f_flgout)); 
#else   
    SetFileOutput(new GFMFileOutput<VectorType,FixupType,FlagType,DIM>(*this)); 
#endif
    SetFixup(new FixupSpecific(integ));
    SetFlagging(new FlaggingSpecific(*this)); 
    AddGFM(new GhostFluidMethod<VectorType,DIM>(
              new F77GFMBoundary<VectorType,DIM>(f_ibndrfl,f_itrans),
              new F77GFMLevelSet<DataType,DIM>(f_lset)));
    _Interpolation = new interpolation_type(*this);
  }  
 
  ~SolverSpecific() {
    DeleteGFM(_GFM[0]);
    delete _Flagging;
    delete _Fixup;
    delete _Interpolation;
  }

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

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

  // 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 % OutputDragEvery != 0) return cfl;

    int me = MY_PROC; 
    double pi = 4.0*std::atan(1.0);
    
    // Read sphere info from F77 common block
    int NSph, Np[5];
    double xm[10], rd[5];
    point_type sum[5];
    f_combl_read(NSph,xm,rd,Np);

    // 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]);
    }

    for (int ns=0; ns<NSph; ns++) {
      int Npoints = 4*Np[ns];
      DataType* p = new DataType[Npoints];
      point_type* xc = new point_type[Npoints];
      point_type* normals = new point_type[Npoints];
      DataType* dist = (DataType*) 0;    
      double dalpha = 0.5*pi/Np[ns];

      register int n;

      // Calculate equidistant points on shere's surfaces
      for (n=0; n<Npoints; n++) {
        double alpha = dalpha*(n+0.5);
        xc[n](0) = xm[base::Dim()*ns  ]+rd[ns]*std::cos(alpha); 
        xc[n](1) = xm[base::Dim()*ns+1]+rd[ns]*std::sin(alpha);
      }
    
      // Interpolate/extrapolate pressure values and assemble results on processor VizServer
      _Interpolation->PointsValues(base::Work(),base::t[0],Npoints,xc,p);
      _Interpolation->ArrayCombine(VizServer,Npoints,p);

      // Approximate normals and assemble results on processor VizServer
      _Interpolation->PointsGeom(base::GFM(0).Phi(),base::t[0],Npoints,xc,normals,dist);
      _Interpolation->ArrayCombine(VizServer,base::Dim()*Npoints,normals[0].data());

      // Do the integration on the surface only on processor VizServer
      if (me == VizServer) {
        sum[ns] = 0.0;
        double ds = 2.*pi*rd[ns]/Npoints; 
        for (n=0; n<Npoints; n++)
          sum[ns] += p[n]*normals[n]*ds; 
      }

      delete [] p;
      delete [] xc;
      delete [] normals;
    }

    // Append to output file only on processor VizServer
    if (me == VizServer) {
      std::ofstream outfile;
      std::ostream* out;
      std::string fname = "drag.dat";
      outfile.open(fname.c_str(), std::ios::out | std::ios::app);
      out = new std::ostream(outfile.rdbuf());
      *out << base::t[0];
      for (int ns=0; ns<NSph; ns++)
        for (int d=0; d<base::Dim(); d++)
          *out << " " << sum[ns](d); 
      *out << std::endl;
      outfile.close();
      delete out;      
    }

    return cfl;
  } 

protected:
  int OutputDragEvery;
  interpolation_type* _Interpolation;