clawpack/applications/euler/2d/SphereLiftOff/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"

#define f_combl_read FORTRAN_NAME(cblread, CBLREAD)
#define f_combl_write FORTRAN_NAME(cblwrite, CBLWRITE)
#define f_ipvel FORTRAN_NAME(ipvel, IPVEL)

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

#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) {
    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,f_ipvel),
              new F77GFMLevelSet<DataType,DIM>(f_lset)));
    _Interpolation = new interpolation_type(*this);
  }  
 
  ~SolverSpecific() {
    DeleteGFM(_GFM[0]);
    delete _Flagging;
    delete _Fixup;
    delete _Interpolation;
  }

  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) {

    int me = MY_PROC; 

    point_type force_old[5];
    double t_old = base::t[0];
    CalculateForce(force_old);

    double cfl = base::Tick(VariableTimeStepping, dtv, cflv, Rejections);

    point_type a, force_new[5];
    double dt = base::t[0]-t_old;
    CalculateForce(force_new);
    
    int NSph, Np[5];
    double xm[10], vm[10], cm[5], rd[5];
    f_combl_read(NSph,xm,vm,cm,rd,Np);

    for (int ns=0; ns<NSph; ns++) {
      a = 0.5/cm[ns]*(force_old[ns]+force_new[ns]);
      vm[base::Dim()*ns  ] += a(0)*dt;
      vm[base::Dim()*ns+1] += a(1)*dt;
      xm[base::Dim()*ns  ] += vm[base::Dim()*ns  ]*dt;
      xm[base::Dim()*ns+1] += vm[base::Dim()*ns+1]*dt;
    }

    f_combl_write(NSph,xm,vm,cm,rd,Np);
    
    // 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++) {
        int d;
        for (d=0; d<base::Dim(); d++)
          *out << " " << force_new[ns](d); 
        for (d=0; d<base::Dim(); d++)
          *out << " " << xm[base::Dim()*ns+d]; 
        for (d=0; d<base::Dim(); d++)
          *out << " " << vm[base::Dim()*ns+d]; 
      }
      *out << std::endl;
      outfile.close();
      delete out;      
    }

    return cfl;
  } 


  void CalculateForce(point_type* sum) {
    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], vm[10], cm[5], rd[5];
    f_combl_read(NSph,xm,vm,cm,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) {
        int Nerr = 0;
        for (n=0; n<Npoints; n++)
          if (p[n]<=-1.e37) Nerr++;
        Npoints -= Nerr;
        sum[ns] = 0.0;
        double ds = 2.*pi*rd[ns]/Npoints; 
        for (n=0; n<Npoints; n++)
          if (p[n]>-1.e37)
            sum[ns] += p[n]*normals[n]*ds; 
      }

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

#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)*5;
      int tag = 10001;
      if (me == VizServer) 
        for (int proc=0; proc<num; proc++) {
          if (proc==VizServer) continue;
          R = MPI_Send((void *)sum, size, MPI_BYTE, proc, tag, comm_service::comm());
          if ( MPI_SUCCESS != R ) 
            comm_service::error_die( "SumCombine", "MPI_Send", R );
        }
      else {
        R = MPI_Recv((void *)sum, size, MPI_BYTE, VizServer, tag, comm_service::comm(), &status);
        if ( MPI_SUCCESS != R ) 
          comm_service::error_die( "ArrayCombine", "MPI_Recv", R );
      }
    }
#endif
  }    

protected:
  interpolation_type* _Interpolation;