vtf-logo

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

00001 // -*- C++ -*-
00002 
00003 // Copyright (C) 2003-2007 California Institute of Technology
00004 // Ralf Deiterding, ralf@amroc.net
00005 
00006 #ifndef AMROC_PROBLEM_H
00007 #define AMROC_PROBLEM_H
00008 
00009 #include <cmath>
00010 #include "euler3.h"
00011 #include "ClpProblem.h"
00012 
00013 #define f_combl_read FORTRAN_NAME(cblread, CBLREAD)
00014 #define f_combl_write FORTRAN_NAME(cblwrite, CBLWRITE)
00015 #define f_ipvel FORTRAN_NAME(ipvel, IPVEL)
00016 
00017 extern "C" {
00018   void f_combl_read(INTEGER& NSph, DOUBLE xm[], DOUBLE vm[], DOUBLE cm[], DOUBLE rd[], INTEGER Np[]); 
00019   void f_combl_write(INTEGER& NSph, DOUBLE xm[], DOUBLE vm[], DOUBLE cm[], DOUBLE rd[], INTEGER Np[]); 
00020   void f_ipvel(); 
00021 }
00022 
00023 #define OWN_GFMAMRSOLVER
00024 #include "ClpStdGFMProblem.h"
00025 #include "AMRGFMInterpolation.h"
00026 
00027 class SolverSpecific : 
00028   public AMRGFMSolver<VectorType,FixupType,FlagType,DIM> {
00029   typedef VectorType::InternalDataType DataType;
00030   typedef AMRGFMSolver<VectorType,FixupType,FlagType,DIM> base;
00031   typedef AMRGFMInterpolation<VectorType,FixupType,FlagType,DIM> interpolation_type;
00032   typedef interpolation_type::point_type point_type;
00033   typedef F77GFMFileOutput<VectorType,FixupType,FlagType,DIM> output_type;
00034 public:
00035   SolverSpecific(IntegratorSpecific& integ, 
00036                  base::initial_condition_type& init,
00037                  base::boundary_conditions_type& bc) : base(integ, init, bc) {
00038     SetLevelTransfer(new F77LevelTransfer<VectorType,DIM>(f_prolong, f_restrict));
00039 #ifdef f_flgout
00040     SetFileOutput(new F77GFMFileOutput<VectorType,FixupType,FlagType,DIM>(*this,f_flgout)); 
00041 #else   
00042     SetFileOutput(new GFMFileOutput<VectorType,FixupType,FlagType,DIM>(*this)); 
00043 #endif
00044     SetFixup(new FixupSpecific(integ));
00045     SetFlagging(new FlaggingSpecific(*this)); 
00046     AddGFM(new GhostFluidMethod<VectorType,DIM>(
00047               new F77GFMBoundary<VectorType,DIM>(f_ibndrfl,f_itrans,f_ipvel),
00048               new F77GFMLevelSet<DataType,DIM>(f_lset)));
00049     _Interpolation = new interpolation_type(*this);
00050   }  
00051  
00052   ~SolverSpecific() {
00053     DeleteGFM(_GFM[0]);
00054     delete _Flagging;
00055     delete _Fixup;
00056     delete _Interpolation;
00057   }
00058 
00059   virtual void SetupData() {
00060     base::SetupData();
00061     _Interpolation->SetupData(base::PGH(), base::NGhosts());
00062   }
00063 
00064   // Overloading of Tick() to evaluate force integral on spheres after each
00065   // level 0 time step.
00066 
00067   virtual double Tick(int VariableTimeStepping, const double dtv[], 
00068                       const double cflv[], int& Rejections) {
00069 
00070     int me = MY_PROC; 
00071 
00072     point_type force_old[5];
00073     double t_old = base::t[0];
00074     CalculateForce(force_old);
00075 
00076     double cfl = base::Tick(VariableTimeStepping, dtv, cflv, Rejections);
00077 
00078     point_type a, force_new[5];
00079     double dt = base::t[0]-t_old;
00080     CalculateForce(force_new);
00081     
00082     int NSph, Np[5];
00083     double xm[15], vm[15], cm[5], rd[5];
00084     f_combl_read(NSph,xm,vm,cm,rd,Np);
00085 
00086     for (int ns=0; ns<NSph; ns++) {
00087       a = 0.5/cm[ns]*(force_old[ns]+force_new[ns]);
00088       vm[base::Dim()*ns  ] += a(0)*dt;
00089       vm[base::Dim()*ns+1] += a(1)*dt;
00090       vm[base::Dim()*ns+2] += a(2)*dt;
00091       xm[base::Dim()*ns  ] += vm[base::Dim()*ns  ]*dt;
00092       xm[base::Dim()*ns+1] += vm[base::Dim()*ns+1]*dt;
00093       xm[base::Dim()*ns+2] += vm[base::Dim()*ns+2]*dt;
00094     }
00095 
00096     f_combl_write(NSph,xm,vm,cm,rd,Np);
00097     
00098     // Append to output file only on processor VizServer
00099     if (me == VizServer) {
00100       std::ofstream outfile;
00101       std::ostream* out;
00102       std::string fname = "drag.dat";
00103       outfile.open(fname.c_str(), std::ios::out | std::ios::app);
00104       out = new std::ostream(outfile.rdbuf());
00105       *out << base::t[0];
00106       for (int ns=0; ns<NSph; ns++) {
00107         int d;
00108         for (d=0; d<base::Dim(); d++)
00109           *out << " " << force_new[ns](d); 
00110         for (d=0; d<base::Dim(); d++)
00111           *out << " " << xm[base::Dim()*ns+d]; 
00112         for (d=0; d<base::Dim(); d++)
00113           *out << " " << vm[base::Dim()*ns+d]; 
00114       }
00115       *out << std::endl;
00116       outfile.close();
00117       delete out;      
00118     }
00119 
00120     return cfl;
00121   } 
00122 
00123 
00124   void CalculateForce(point_type* sum) {
00125     int me = MY_PROC; 
00126 
00127     double pi = 4.0*std::atan(1.0);
00128     
00129     // Read sphere info from F77 common block
00130     int NSph, Np[5];
00131     double xm[15], vm[15], cm[5], rd[5];
00132     f_combl_read(NSph,xm,vm,cm,rd,Np);
00133 
00134     // Calculate pressure
00135     for (register int l=0; l<=FineLevel(base::GH()); l++) {
00136       int Time = CurrentTime(base::GH(),l);
00137       int press_cnt = base::Dim()+4;      
00138       ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time, l, 
00139                                               press_cnt, base::t[l]);
00140     }
00141 
00142     for (int ns=0; ns<NSph; ns++) {
00143       int Npoints = 8*Np[ns]*Np[ns];
00144       DataType* p = new DataType[Npoints];
00145       point_type* xc = new point_type[Npoints];
00146       point_type* normals = new point_type[Npoints];
00147       DataType* ds = new DataType[Npoints];
00148       DataType* dist = (DataType*) 0;    
00149       double dalpha = 0.5*pi/Np[ns];
00150       double dbeta = 0.5*pi/Np[ns];
00151       register int n = 0;
00152 
00153       // Calculate points on shere's surfaces
00154       for (register int l=0; l<4*Np[ns]; l++) {
00155         double alpha = dalpha*(l+0.5);
00156         for (register int m=0; m<2*Np[ns]; m++) {
00157           double beta = dalpha*(m+0.5);
00158           xc[n](0) = xm[base::Dim()*ns  ]+rd[ns]*std::cos(alpha)*std::sin(beta); 
00159           xc[n](1) = xm[base::Dim()*ns+1]+rd[ns]*std::sin(alpha)*std::sin(beta);
00160           xc[n](2) = xm[base::Dim()*ns+2]+rd[ns]*std::cos(beta);
00161           ds[n] = pi*rd[ns]*rd[ns]*dbeta*std::sin(beta)/(2.*Np[ns]);
00162           n++;
00163         }
00164       }
00165     
00166       // Interpolate/extrapolate pressure values and assemble results on processor VizServer
00167       _Interpolation->PointsValues(base::Work(),base::t[0],Npoints,xc,p);
00168       _Interpolation->ArrayCombine(VizServer,Npoints,p);
00169 
00170       // Approximate normals and assemble results on processor VizServer
00171       _Interpolation->PointsGeom(base::GFM(0).Phi(),base::t[0],Npoints,xc,normals,dist);
00172       _Interpolation->ArrayCombine(VizServer,base::Dim()*Npoints,normals[0].data());
00173 
00174       // Do the integration on the surface only on processor VizServer
00175       if (me == VizServer) {
00176         sum[ns] = 0.0; 
00177         for (n=0; n<Npoints; n++)
00178           if (p[n]>-1.e37) 
00179             sum[ns] += p[n]*normals[n]*ds[n]; 
00180       }
00181 
00182       delete [] p;
00183       delete [] xc;
00184       delete [] normals;
00185       delete [] ds;
00186     }
00187 
00188 #ifdef DAGH_NO_MPI
00189 #else
00190     if (comm_service::dce() && comm_service::proc_num() > 1) {  
00191       int num = comm_service::proc_num();
00192       MPI_Status status;
00193       int R;
00194       int size = sizeof(DataType)*5;
00195       int tag = 10001;
00196       if (me == VizServer) 
00197         for (int proc=0; proc<num; proc++) {
00198           if (proc==VizServer) continue;
00199           R = MPI_Send((void *)sum, size, MPI_BYTE, proc, tag, comm_service::comm());
00200           if ( MPI_SUCCESS != R ) 
00201             comm_service::error_die( "SumCombine", "MPI_Send", R );
00202         }
00203       else {
00204         R = MPI_Recv((void *)sum, size, MPI_BYTE, VizServer, tag, comm_service::comm(), &status);
00205         if ( MPI_SUCCESS != R ) 
00206           comm_service::error_die( "ArrayCombine", "MPI_Recv", R );
00207       }
00208     }
00209 #endif
00210   }    
00211 
00212 protected:
00213   interpolation_type* _Interpolation;