clawpack/applications/euler/3d/MovingSpheres/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 "euler3.h"
#include "ClpProblem.h"

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

extern "C" {
  void f_combl_write(INTEGER& NSph, DOUBLE* xm, DOUBLE* vm, DOUBLE rd[]); 
  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), 
    CouplingLevel(0), NSph(0), Np(0), xm(0), vm(0), fm(0), cm(0), rd(0), tm(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)));
    _Interpolation = new interpolation_type(*this);
    pi = 4.0*std::atan(1.0);
  }  
 
  ~SolverSpecific() {
    DeleteGFM(_GFM[0]);
    delete _Flagging;
    delete _Fixup;
    delete _Interpolation;
    if (Np) delete [] Np;
    if (xm) delete [] xm;
    if (vm) delete [] vm;
    if (fm) delete [] fm;
    if (cm) delete [] cm;
    if (rd) delete [] rd;
    if (tm) delete [] tm;
  }

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

  virtual void SetupData() {
    base::SetupData();
    _Interpolation->SetupData(base::PGH(), base::NGhosts());
    if (CouplingLevel>base::MaxLev-1) 
      CouplingLevel=base::MaxLev-1;    
  }

  virtual void Initialize(double& t, double &dt) {
    std::ifstream ifs("spheres.dat", std::ios::in);
    ifs >> NSph;
    Np = new int[NSph];
    xm = new point_type[NSph];
    vm = new point_type[NSph];
    fm = new point_type[NSph];
    cm = new DataType[NSph];
    rd = new DataType[NSph];
    tm = new DataType[NSph];
    for (int n=0; n<NSph; n++) {
      ifs >> xm[n](0) >> xm[n](1) >> xm[n](2) >> rd[n] 
          >> vm[n](0) >> vm[n](1) >> vm[n](2) >> cm[n]
          >> tm[n] >> Np[n];
      cm[n] = cm[n]*pi*4./3.*rd[n]*rd[n]*rd[n];
    }
    f_combl_write(NSph,&xm[0](0),&vm[0](0),rd);
    ifs.close();
    base::Initialize(t,dt);
    OutputFiles(t);
  }

  virtual void Restart(double& t, double& dt) {
    std::ifstream ifs("spheres.cp", std::ios::in);
    ifs.read((char*)&NSph,sizeof(int));
    Np = new int[NSph];
    xm = new point_type[NSph];
    vm = new point_type[NSph];
    fm = new point_type[NSph];
    cm = new DataType[NSph];
    rd = new DataType[NSph];
    tm = new DataType[NSph];
    ifs.read((char*)Np,NSph*sizeof(int));    
    ifs.read((char*)xm,NSph*sizeof(point_type));
    ifs.read((char*)vm,NSph*sizeof(point_type));
    ifs.read((char*)fm,NSph*sizeof(point_type));
    ifs.read((char*)cm,NSph*sizeof(DataType));
    ifs.read((char*)rd,NSph*sizeof(DataType));
    ifs.read((char*)tm,NSph*sizeof(DataType));
    ifs.close();
    f_combl_write(NSph,&xm[0](0),&vm[0](0),rd);
    base::Restart(t,dt);
  }
  virtual void Checkpointing() {
    int me = MY_PROC; 
    if (me == VizServer) {
      std::ofstream ofs("spheres.cp", std::ios::out);
      ofs.write((char*)&NSph,sizeof(int));
      ofs.write((char*)Np,NSph*sizeof(int));    
      ofs.write((char*)xm,NSph*sizeof(point_type));
      ofs.write((char*)vm,NSph*sizeof(point_type));
      ofs.write((char*)fm,NSph*sizeof(point_type));
      ofs.write((char*)cm,NSph*sizeof(DataType));
      ofs.write((char*)rd,NSph*sizeof(DataType));
      ofs.write((char*)tm,NSph*sizeof(DataType));
      ofs.close();
    }
    base::Checkpointing();
  }

  virtual double IntegrateLevel(vec_grid_fct_type& u, const int Time, const int Level, 
                                double t, double dt, bool DoFixup, double tc, 
                                const int which) {

    double cfl = base::IntegrateLevel(u,Time,Level,t,dt,DoFixup,tc,which);
#ifdef DEBUG_PRINT
    ( comm_service::log() << "*** base::IntegrateLevel() " << Level << " finished.\n" ).flush();
#endif

    if (which==0 && Level==CouplingLevel) {
      point_type* fm_old = new point_type[NSph];
      for (int ns=0; ns<NSph; ns++) 
        fm_old[ns] = fm[ns];
      point_type a;
      CalculateForce();   
      for (int ns=0; ns<NSph; ns++) {
        if (t+dt<=tm[ns]) continue;
        a = 0.5/cm[ns]*(fm[ns]+fm_old[ns]);
        vm[ns] += a*dt;
        xm[ns] += vm[ns]*dt;
      }
      delete [] fm_old;
      f_combl_write(NSph,&xm[0](0),&vm[0](0),rd);
      OutputFiles(t+dt);
    }

    return cfl;
  }  

  void CalculateForce() {
    // 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]);
    }

    double *force = new double[3*NSph];
    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 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[ns](0)+rd[ns]*std::cos(alpha)*std::sin(beta); 
          xc[n](1) = xm[ns](1)+rd[ns]*std::sin(alpha)*std::sin(beta);
          xc[n](2) = xm[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 normals
      _Interpolation->PointsValues(base::Work(),base::t[0],Npoints,xc,p);
      _Interpolation->PointsGeom(base::GFM(0).Phi(),base::t[0],Npoints,xc,normals,dist);

      force[3*ns  ] = 0.0; 
      force[3*ns+1] = 0.0; 
      force[3*ns+2] = 0.0; 
      for (n=0; n<Npoints; n++)
        if (p[n]>-1.e37) {
          force[3*ns  ] += p[n]*normals[n](0)*ds[n]; 
          force[3*ns+1] += p[n]*normals[n](1)*ds[n]; 
          force[3*ns+2] += p[n]*normals[n](2)*ds[n]; 
        }

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

#ifdef DAGH_NO_MPI
    for (int ns=0; ns<NSph; ns++) {
      fm[ns](0) = force[3*ns  ];
      fm[ns](1) = force[3*ns+1];
      fm[ns](2) = force[3*ns+2];
    }
#else
    int Nval = 3*NSph;
    double *global_force = new double[3*NSph];
    MPI_Allreduce(force, global_force, Nval, MPI_DOUBLE, MPI_SUM, base::Comm());
    for (int ns=0; ns<NSph; ns++) {
      fm[ns](0) = global_force[3*ns  ];
      fm[ns](1) = global_force[3*ns+1];
      fm[ns](2) = global_force[3*ns+2];
    }
#endif

    delete [] force;
  }    

  void OutputFiles(double t) {
    int me = MY_PROC; 
    if (me == VizServer) {
      std::ofstream ofs1("positions.txt", std::ios::out | std::ios::app );
      ofs1 << t;
      for (int ns=0; ns<NSph; ns++)
        ofs1 << " " << xm[ns](0) << " " << xm[ns](1) << " " << xm[ns](2);  
      ofs1 << std::endl;
        ofs1.close();

      std::ofstream ofs2("forces.txt", std::ios::out | std::ios::app );
      ofs2 << t;
      for (int ns=0; ns<NSph; ns++)
        ofs2 << " " << fm[ns](0) << " " << fm[ns](1) << " " << fm[ns](2);  
      ofs2 << std::endl;
        ofs2.close();

      std::ofstream ofs3("velocities.txt", std::ios::out | std::ios::app );
      ofs3 << t;
      for (int ns=0; ns<NSph; ns++)
        ofs3 << " " << vm[ns](0) << " " << vm[ns](1) << " " << vm[ns](2);  
      ofs3 << std::endl;
        ofs3.close();
    }
  }    

protected:
  interpolation_type* _Interpolation;
  int CouplingLevel, NSph, *Np;
  point_type *xm, *vm, *fm;
  DataType *cm, *rd, *tm, pi;