clawpack/applications/eulerm/2d/WaterShockTube/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 "eulerm2.h"

#undef  NAUX
#define NAUX       1

#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(DOUBLE xm[], DOUBLE vm[], DOUBLE& cm, DOUBLE& rd, INTEGER& Np); 
  void f_combl_write(DOUBLE xm[], DOUBLE vm[], DOUBLE& cm, DOUBLE& rd, INTEGER& Np); 
  void f_ipvel(); 
}

#define OUTPUTPOINTSMAX  20
#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;
  typedef GhostFluidMethod<VectorType,DIM>::boundary_conditions_type gfm_bc_type;
public:
  typedef base::vec_grid_data_type vec_grid_data_type;  

  SolverSpecific(IntegratorSpecific& integ, 
                 base::initial_condition_type& init,
                 base::boundary_conditions_type& bc) : base(integ, init, bc), OutputEvery(0), Nxc(0) {
    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),
              (gfm_bc_type*)(new F77BoundaryConditions<Vector<DataType,1>,DIM>(f_boundary))));
    _Interpolation = new interpolation_type(*this);
  }  
 
  ~SolverSpecific() {
    delete _GFM[0]->GetBoundaryConditionsP();
    DeleteGFM(_GFM[0]);
    delete _Flagging;
    delete _Fixup;
    delete _Interpolation;
  }

  virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
    base::register_at(Ctrl,prefix);
    SensorsCtrl = base::LocCtrl.getSubDevice("SensorsOutput");
    RegisterAt(SensorsCtrl,"OutputEvery",OutputEvery);
    RegisterAt(SensorsCtrl,"NSensors",Nxc);
    char VariableName[32];
    for (register int i=0; i<OUTPUTPOINTSMAX; i++) 
      for (register int d=0; d<base::Dim(); d++) {
        std::sprintf(VariableName,"xc(%d,%d)",i+1,d+1);
        RegisterAt(SensorsCtrl,VariableName,xc[i](d));
      }
  } 
  virtual void register_at(ControlDevice& Ctrl) {
    base::register_at(Ctrl);
  }

  virtual void SetupData() {
    base::SetupData();
    _GFM[0]->BoundaryConditions_().SetupData(base::PGH(), base::NGhosts());
    _Interpolation->SetupData(base::PGH(), base::NGhosts());
    if (OutputEvery<0) OutputEvery=0;
    if (Nxc<0) Nxc=0;
    if (Nxc>OUTPUTPOINTSMAX) Nxc=OUTPUTPOINTSMAX;
  }

  // Overloading of Advance() to add pointwise output after each level 0 time step
  virtual void Advance(double& t, double& dt) {

    int me = MY_PROC; 
    point_type force_old;
    double t_old = t;
    CalculateForce(force_old);

    base::Advance(t,dt);

    point_type a, force_new;
    CalculateForce(force_new);
    
    int Np;
    double xm[2], vm[2], cm, rd;
    f_combl_read(xm,vm,cm,rd,Np);

    a = 0.5/cm*(force_old+force_new);
    vm[0] += a(0)*dt;
    vm[1] += a(1)*dt;
    xm[0] += vm[0]*dt;
    xm[1] += vm[1]*dt;

    f_combl_write(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 = "boundary.txt";
      outfile.open(fname.c_str(), std::ios::out | std::ios::app);
      out = new std::ostream(outfile.rdbuf());
      *out << t;
      int d;
      for (d=0; d<base::Dim(); d++)
        *out << " " << force_new(d); 
      for (d=0; d<base::Dim(); d++)
        *out << " " << xm[d]; 
      for (d=0; d<base::Dim(); d++)
        *out << " " << vm[d]; 
      *out << std::endl;
      outfile.close();
      delete out;      
    }

    if (!OutputEvery) return;

    int TimeZero  = CurrentTime(base::GH(),0);
    TimeZero /= RefinedBy(base::GH(),MaxLevel(base::GH()));
    if (TimeZero % OutputEvery != 0) return;

    VectorType* dat = new VectorType[Nxc];

    // Interpolate/extrapolate pressure values and assemble results on processor VizServer
    _Interpolation->PointsValues(base::U(),t,Nxc,xc,dat);
    int NValues = Nxc*base::NEquations();
    _Interpolation->ArrayCombine(VizServer,NValues,&(dat[0](0)));

    // Append to output file only on processor VizServer
    if (MY_PROC == VizServer) {
      DataType* output_data = new DataType[_FileOutput->Ncnt()*Nxc];
      BBox bb(2,0,0,Nxc-1,0,1,1);
      vec_grid_data_type val(bb,dat);
      for (int cnt=1;cnt<=_FileOutput->Ncnt(); cnt++) {
        if (_FileOutput->OutputName(cnt-1).c_str()[0] == '-') continue;
        grid_data_type output(bb,output_data+(cnt-1)*Nxc);
        ((output_type::out_2_func_type) (((output_type*) base::_FileOutput)->GetOutFunc()))
          (FA(2,val),FA(2,output),BOUNDING_BOX(bb),base::NEquations(),cnt,t);
        DataType* dummy;
        output.deallocate(dummy);
      }
      VectorType* dummy;
      val.deallocate(dummy);
      std::ofstream outfile;
      std::ostream* out;
      std::string fname = "sensors.txt";
      outfile.open(fname.c_str(), std::ios::out | std::ios::app);
      out = new std::ostream(outfile.rdbuf());
      *out << t << " ";
      for (register int j=0;j<Nxc; j++)
        for (int cnt=0;cnt<_FileOutput->Ncnt(); cnt++) 
          if (_FileOutput->OutputName(cnt).c_str()[0] != '-') 
            *out << " " << output_data[cnt*Nxc+j];
      *out << std::endl;
      delete [] output_data;
      outfile.close();
      delete out;      
    }
    delete [] dat;
  } 

  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 Np;
    double xm[2], vm[2], cm, rd;
    f_combl_read(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]);
    }

    DataType* p = new DataType[Np];
    point_type* xc = new point_type[Np];

    register int n;

    for (n=0; n<Np; n++) {
      xc[n](0) = xm[0]; 
      xc[n](1) = xm[1]+(2.*n/Np-1.)*rd;
    }
    
    // Interpolate/extrapolate pressure values and assemble results on processor VizServer
    _Interpolation->PointsValues(base::Work(),base::t[0],Np,xc,p);
    _Interpolation->ArrayCombine(VizServer,Np,p);

    // Do the integration on the surface only on processor VizServer
    if (me == VizServer) {
      sum = 0.;
      int c=0;
      for (n=0; n<Np; n++)
        if (p[n]>-1.e37) {
          sum(0) += p[n]; c++;
        }
      sum(0) /= c;
      sum(0) -= 101325.;
      sum(0) *= pi*rd*rd;

      delete [] p;
      delete [] xc;
    }

#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);
      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:
  int OutputEvery, Nxc;
  point_type xc[OUTPUTPOINTSMAX];
  interpolation_type* _Interpolation;