clawpack/applications/euler/3d/Spheres/src/Problem.h

// -*- C++ -*-

// Copyright (C) 2003-2007 California Institute of Technology
// Ralf Deiterding, ralf@amroc.net
// Stuart Laurence

#ifndef AMROC_PROBLEM_H
#define AMROC_PROBLEM_H

#include <cmath>
#include "euler3.h"
#include "ClpProblem.h"

#define f_combl_read FORTRAN_NAME(cblread, CBLREAD)

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

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

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());
    if (OutputDragEvery<=0) OutputDragEvery=1;
  }

  // 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[15], 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 = 8*Np[ns]*Np[ns];
      DataType* p = new DataType[Npoints];
      point_type* xc = new point_type[Npoints];
      point_type* normals = new point_type[Npoints];
      DataType* ds = new DataType[Npoints];
      DataType* dist = (DataType*) 0;    
      double dalpha = 0.5*pi/Np[ns];
      double dbeta = 0.5*pi/Np[ns];
      register int n = 0;

      // Calculate equidistant points on shere's surfaces
      for (register int l=0; l<4*Np[ns]; l++) {
        double alpha = dalpha*(l+0.5);
        for (register int m=0; m<2*Np[ns]; m++) {
          double beta = dalpha*(m+0.5);
          xc[n](0) = xm[base::Dim()*ns  ]+rd[ns]*std::cos(alpha)*std::sin(beta); 
          xc[n](1) = xm[base::Dim()*ns+1]+rd[ns]*std::sin(alpha)*std::sin(beta);
          xc[n](2) = xm[base::Dim()*ns+2]+rd[ns]*std::cos(beta);
          ds[n] = pi*rd[ns]*rd[ns]*dbeta*std::sin(beta)/(2.*Np[ns]);
          n++;
        }
      }
    
      // 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; 
        for (n=0; n<Npoints; n++)
          sum[ns] += p[n]*normals[n]*ds[n]; 
      }

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

    // 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;