fsi/sfc-amroc/WaterBlastPlastic/src/FluidProblem.h

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 
00003 // Copyright (C) 2003-2007 California Institute of Technology
00004 // Ralf Deiterding, ralf@cacr.caltech.edu
00005 
00006 #ifndef AMROC_FLUIDPROBLEM_H
00007 #define AMROC_FLUIDPROBLEM_H
00008 
00009 #include "eulerm3.h"
00010 
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_piston FORTRAN_NAME(ipvelpiston, IPVELPISTON)
00016 #define f_lset_piston FORTRAN_NAME(lspiston, LSPISTON)
00017 
00018 extern "C" {
00019   void f_combl_read(DOUBLE& xm, DOUBLE& vm, DOUBLE& cm, DOUBLE& rd, INTEGER& Np); 
00020   void f_combl_write(DOUBLE& xm, DOUBLE& vm, DOUBLE& cm, DOUBLE& rd, INTEGER& Np); 
00021   void f_ipvel_piston(); 
00022   void f_lset_piston(); 
00023 }
00024 
00025 #define MaxDPoints (10)
00026 #define OUTPUTPOINTSMAX  20
00027 
00028 #define OWN_FLAGGING
00029 #define OWN_ELCGFMAMRSOLVER
00030 #include "ClpStdELCGFMProblem.h" 
00031 #include "AMRGFMInterpolation.h"
00032 
00033 
00034 class FlaggingSpecific : 
00035   public AMRFlagging<VectorType,FixupType,FlagType,DIM> {
00036   typedef AMRFlagging<VectorType,FixupType,FlagType,DIM> base;
00037 public:
00038   FlaggingSpecific(base::solver_type& solver) : base(solver) {
00039       base::AddCriterion(new ByValue<VectorType,FlagType,DIM>());
00040       base::AddCriterion(new ScaledGradient<VectorType,FlagType,DIM>());
00041       base::AddCriterion(new LimiterType<VectorType,FlagType,DIM>());
00042       base::AddCriterion(new AbsoluteError<VectorType,FixupType,FlagType,DIM>(solver));
00043       base::AddCriterion(new RelativeError<VectorType,FixupType,FlagType,DIM>(solver));
00044 #ifdef f_flgout
00045       base::AddCriterion(new F77ByValue<VectorType,FlagType,DIM>(f_flgout));
00046       base::AddCriterion(new F77ScaledGradient<VectorType,FlagType,DIM>(f_flgout));
00047       base::AddCriterion(new F77LimiterType<VectorType,FlagType,DIM>(f_flgout));
00048       base::AddCriterion(new F77AbsoluteError<VectorType,FixupType,FlagType,DIM>(solver,f_flgout));
00049       base::AddCriterion(new F77RelativeError<VectorType,FixupType,FlagType,DIM>(solver,f_flgout));
00050 #endif
00051       dcmin = new int[MaxDPoints*base::Dim()];
00052       dcmax = new int[MaxDPoints*base::Dim()];
00053       dclev = new int[MaxDPoints];
00054       for (int d=0; d<MaxDPoints*base::Dim(); d++) {
00055         dcmin[d] = 0; dcmax[d] = -1;
00056       }
00057       for (int i=0; i<MaxDPoints; i++)
00058         dclev[i] = 1000000;
00059     }
00060 
00061   ~FlaggingSpecific() { 
00062     DeleteAllCriterions(); 
00063     delete [] dcmin;
00064     delete [] dcmax;
00065   }
00066 
00067   virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00068     base::register_at(Ctrl, prefix);
00069     char VariableName[32];
00070     for (int i=0; i<MaxDPoints; i++) {
00071       for (int d=0; d<base::Dim(); d++) {
00072         sprintf(VariableName,"DCellsMin(%d,%d)",i+1,d+1);
00073         RegisterAt(base::LocCtrl,VariableName,dcmin[i*base::Dim()+d]);
00074         sprintf(VariableName,"DCellsMax(%d,%d)",i+1,d+1);
00075         RegisterAt(base::LocCtrl,VariableName,dcmax[i*base::Dim()+d]);
00076       }
00077       sprintf(VariableName,"DMinLevel(%d)",i+1);
00078       RegisterAt(base::LocCtrl,VariableName,dclev[i]);
00079     }
00080   }
00081   virtual void register_at(ControlDevice& Ctrl) {
00082     register_at(Ctrl, "");
00083   }
00084 
00085   virtual void SetFlags(const int Time, const int Level, double t, double dt) {
00086     base::SetFlags(Time, Level, t, dt);
00087     for (int i=0; i<MaxDPoints; i++) 
00088       if (Level>=dclev[i]) {
00089         Coords lb(base::Dim(), dcmin+i*base::Dim()), ub(base::Dim(), dcmax+i*base::Dim()), 
00090                   s(base::Dim(),1);
00091         BBox bbox(lb*RefinedBy(base::GH(),MaxLevel(base::GH())),
00092                   ub*RefinedBy(base::GH(),MaxLevel(base::GH())),
00093                   s*StepSize(base::GH(),Level));
00094         forall (base::Flags(),Time,Level,c)  
00095           base::Flags()(Time,Level,c).equals(GoodPoint, bbox);
00096         end_forall    
00097       }
00098   }
00099 
00100 private:
00101   int *dcmin, *dcmax, *dclev;
00102 };
00103 
00104 
00105 template <class DataType, int dim>
00106 class F77CPTLevelSet : public CPTLevelSet<DataType,dim> {
00107   typedef CPTLevelSet<DataType,dim> base;
00108 public:
00109   typedef typename base::grid_fct_type grid_fct_type;  
00110   typedef typename base::grid_data_type grid_data_type;
00111   typedef generic_fortran_func generic_func_type; 
00112 
00113   typedef void (*lset_3_func_type) ( const INTEGER& maxmx, const INTEGER& maxmy, const INTEGER& maxmz,
00114                                      const INTEGER& mbc,
00115                                      const INTEGER& mx, const INTEGER& my, const INTEGER& mz,
00116                                      const DOUBLE x[], const DOUBLE y[], const DOUBLE z[],
00117                                      const DOUBLE& dx, const DOUBLE& dy, const DOUBLE& dz, 
00118                                      DataType q[], const DOUBLE& t);
00119 
00120   F77CPTLevelSet() : base(), f_lset(0) {}
00121   F77CPTLevelSet(generic_func_type lset) : base(), f_lset(lset) {}
00122 
00123   virtual void SetPhi(grid_fct_type& phi, const int Time, const int Level, double t) { 
00124     base::SetPhi(phi,Time,Level,t);
00125     forall (phi,Time,Level,c)
00126       SetGrid(phi(Time,Level,c),Level,t);    
00127     end_forall
00128   }
00129 
00130   virtual void SetGrid(grid_data_type& gdphi, const int& level, const double& t) {   
00131     assert(f_lset != (generic_func_type) 0);
00132     Coords ex = gdphi.extents();
00133     DCoords lbcorner = base::GH().worldCoords(gdphi.lower(), gdphi.stepsize());
00134     DCoords dx = base::GH().worldStep(gdphi.stepsize());
00135     int maxm[3], mx[3], d;
00136     DataType* x[3];
00137     for (d=0; d<dim; d++) {
00138       maxm[d] = ex(d);
00139       mx[d] = ex(d)-2*base::NGhosts();
00140       x[d] = new DataType[maxm[d]];
00141       for (int i=0; i<maxm[d]; i++) x[d][i] = (i+0.5)*dx(d)+lbcorner(d);
00142     }
00143 
00144     ((lset_3_func_type) f_lset)(AA(3,mx), base::NGhosts(), AA(3,mx), 
00145                                 AA(3,x), AA(3,dx), gdphi.data(), t);
00146 
00147     for (d=0; d<dim; d++) 
00148       delete [] x[d];
00149   }
00150 
00151   inline void SetFunc(generic_func_type lset) { f_lset = lset; }
00152   generic_func_type GetFunc() const { return f_lset; }
00153 
00154 protected:
00155   generic_func_type f_lset;
00156 };
00157 
00158 
00159 class FluidSolverSpecific : 
00160   public AMRELCGFMSolver<VectorType,FixupType,FlagType,DIM> {
00161   typedef VectorType::InternalDataType DataType;
00162   typedef AMRELCGFMSolver<VectorType,FixupType,FlagType,DIM> base;
00163   typedef AMRGFMInterpolation<VectorType,FixupType,FlagType,DIM> interpolation_type;
00164 public:
00165   typedef F77GFMFileOutput<VectorType,FixupType,FlagType,DIM> output_type;
00166   typedef base::point_type point_type;
00167   typedef base::eul_comm_type eul_comm_type;
00168 
00169   FluidSolverSpecific() : base(_IntegratorSpecific, _InitialConditionSpecific, 
00170                                _BoundaryConditionsSpecific) {
00171     SetLevelTransfer(new F77LevelTransfer<VectorType,DIM>(f_prolong, f_restrict));
00172     SetFileOutput(new output_type(*this,f_flgout)); 
00173     SetFixup(new FixupSpecific(_IntegratorSpecific));
00174     SetFlagging(new FlaggingSpecific(*this)); 
00175     AddGFM(new GhostFluidMethod<VectorType,DIM>(
00176               new F77ELCGFMBoundary<VectorType,DIM>(f_ibndrfl,f_itrans,*this),
00177               new F77CPTLevelSet<DataType,DIM>(f_lset)));
00178     AddGFM(new GhostFluidMethod<VectorType,DIM>(
00179               new F77GFMBoundary<VectorType,DIM>(f_ibndrfl,f_itrans,f_ipvel_piston),
00180               new F77GFMLevelSet<DataType,DIM>(f_lset_piston)));
00181     SetCoupleGFM(0);
00182     _Interpolation = new interpolation_type(*this);
00183   }  
00184  
00185   ~FluidSolverSpecific() {
00186     DeleteGFM(_GFM[1]);
00187     DeleteGFM(_GFM[0]);
00188     delete _Flagging;
00189     delete _Fixup;
00190     delete _Interpolation;
00191   }
00192 
00193   virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00194     base::register_at(Ctrl,prefix);
00195     SensorsCtrl = base::LocCtrl.getSubDevice("SensorsOutput");
00196     RegisterAt(SensorsCtrl,"OutputEvery",OutputEvery);
00197     RegisterAt(SensorsCtrl,"NSensors",Nxc);
00198     char VariableName[32];
00199     for (register int i=0; i<OUTPUTPOINTSMAX; i++) 
00200       for (register int d=0; d<base::Dim(); d++) {
00201         std::sprintf(VariableName,"xc(%d,%d)",i+1,d+1);
00202         RegisterAt(SensorsCtrl,VariableName,xc[i](d));
00203       }
00204   } 
00205   virtual void register_at(ControlDevice& Ctrl) {
00206     base::register_at(Ctrl);
00207   }
00208 
00209   virtual void SetupData() {
00210     base::SetupData();
00211     _Interpolation->SetupData(base::PGH(), base::NGhosts());
00212     if (OutputEvery<0) OutputEvery=0;
00213     if (Nxc<0) Nxc=0;
00214     if (Nxc>OUTPUTPOINTSMAX) Nxc=OUTPUTPOINTSMAX;
00215   }
00216 
00217   virtual void Restart(double& t, double& dt) {
00218     std::ifstream ifs("boundary.cp", std::ios::in);
00219     int Np;
00220     double xm, vm, cm, rd;
00221     f_combl_read(xm,vm,cm,rd,Np);
00222     ifs >> xm >> vm;
00223     ifs.close();
00224     f_combl_write(xm,vm,cm,rd,Np);
00225     base::Restart(t,dt);
00226   }
00227   virtual void Checkpointing() {
00228     int me = MY_PROC; 
00229     if (me == VizServer) {
00230       int Np;
00231       double xm, vm, cm, rd;
00232       f_combl_read(xm,vm,cm,rd,Np);
00233       std::ofstream ofs("boundary.cp", std::ios::out);
00234       ofs << xm << " " << vm << std::endl; 
00235       ofs.close();
00236     }
00237     base::Checkpointing();
00238   }
00239 
00240   virtual void SendBoundaryData() {
00241     START_WATCH
00242       for (register int l=0; l<=FineLevel(base::GH()); l++) {
00243         int Time = CurrentTime(base::GH(),l);
00244         ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time,
00245                                                 l, base::Dim()+4, base::t[l]);
00246         if (CurrentTime(base::GH(),base::CouplingLevel) != Time) 
00247           ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time+TimeStep(base::U(),l),
00248                                                   l, base::Dim()+4, base::t[l]+base::dt[l]);      
00249       }
00250     END_WATCH(FLUID_CPL_PRESSURE_CALCULATE)
00251     base::SendBoundaryData();
00252   }
00253 
00254   virtual void ModifyELCSendData(eul_comm_type* eulComm) {
00255     DataType distance = -GFM(CoupleGFM()).Boundary().Cutoff();
00256     int Npoints = _eulComm->getNumberOfFaces();
00257     double *press_data = _eulComm->getPressuresData();
00258     const point_type *norm = reinterpret_cast<const point_type *>(_eulComm->getFaceNormalsData());
00259     const point_type *cent = reinterpret_cast<const point_type *>(_eulComm->getFaceCentroidsData());
00260     for (register int n=0; n<Npoints; n++) {
00261       if (press_data[n]!=std::numeric_limits<DataType>::max()) {
00262         press_data[n] *= -1; press_data[n] += base::PressureConversionFactor()*101325.; 
00263       }
00264       point_type p = cent[n];
00265       if (distance>1.e-12) p += distance*norm[n];    
00266       if (std::sqrt(p(1)*p(1)+p(2)*p(2))>0.032)
00267         press_data[n] = std::numeric_limits<DataType>::max();
00268     }
00269   }
00270 
00271   // Overloading of Advance() to add pointwise output after each level 0 time step
00272   // and propagate piston
00273   virtual void Advance(double& t, double& dt) {
00274 
00275     int me = MY_PROC; 
00276     DataType force_old;
00277     if (base::GFM(1).IsUsed())
00278       CalculateForce(force_old);
00279 
00280     base::Advance(t,dt);
00281 
00282     if (base::GFM(1).IsUsed()) {
00283       DataType a, force_new;
00284       CalculateForce(force_new);
00285       
00286       int Np;
00287       double xm, vm, cm, rd;
00288       f_combl_read(xm,vm,cm,rd,Np);
00289       
00290       a = 0.5/cm*(force_old+force_new);
00291       vm += a*dt;
00292       xm += vm*dt;
00293       
00294       f_combl_write(xm,vm,cm,rd,Np);
00295     
00296       // Append to output file only on processor VizServer
00297       if (me == VizServer) {
00298         std::ofstream outfile;
00299         std::ostream* out;
00300         std::string fname = "boundary.txt";
00301         outfile.open(fname.c_str(), std::ios::out | std::ios::app);
00302         out = new std::ostream(outfile.rdbuf());
00303         *out << t << " " << force_new << " " << xm << " " << vm << std::endl;
00304         outfile.close();
00305         delete out;      
00306       }
00307     }
00308 
00309     if (!OutputEvery) return;
00310 
00311     int TimeZero  = CurrentTime(base::GH(),0);
00312     TimeZero /= RefinedBy(base::GH(),MaxLevel(base::GH()));
00313     if (TimeZero % OutputEvery != 0) return;
00314 
00315     VectorType* dat = new VectorType[Nxc];
00316 
00317     // Interpolate/extrapolate pressure values and assemble results on processor VizServer
00318     _Interpolation->PointsValues(base::U(),t,Nxc,xc,dat);
00319     int NValues = Nxc*base::NEquations();
00320     _Interpolation->ArrayCombine(VizServer,NValues,&(dat[0](0)));
00321 
00322     // Append to output file only on processor VizServer
00323     if (MY_PROC == VizServer) {
00324       DataType* output_data = new DataType[_FileOutput->Ncnt()*Nxc];
00325       BBox bb(3,0,0,0,Nxc-1,0,0,1,1,1);
00326       vec_grid_data_type val(bb,dat);
00327       for (int cnt=1;cnt<=_FileOutput->Ncnt(); cnt++) {
00328         if (_FileOutput->OutputName(cnt-1).c_str()[0] == '-') continue;
00329         grid_data_type output(bb,output_data+(cnt-1)*Nxc);
00330         ((output_type::out_3_func_type) (((output_type*) base::_FileOutput)->GetOutFunc()))
00331           (FA(3,val),FA(3,output),BOUNDING_BOX(bb),base::NEquations(),cnt,t);
00332         DataType* dummy;
00333         output.deallocate(dummy);
00334       }
00335       VectorType* dummy;
00336       val.deallocate(dummy);
00337       std::ofstream outfile;
00338       std::ostream* out;
00339       std::string fname = "sensors.txt";
00340       outfile.open(fname.c_str(), std::ios::out | std::ios::app);
00341       out = new std::ostream(outfile.rdbuf());
00342       *out << t << " ";
00343       for (register int j=0;j<Nxc; j++)
00344         for (int cnt=0;cnt<_FileOutput->Ncnt(); cnt++) 
00345           if (_FileOutput->OutputName(cnt).c_str()[0] != '-') 
00346             *out << " " << output_data[cnt*Nxc+j];
00347       *out << std::endl;
00348       delete [] output_data;
00349       outfile.close();
00350       delete out;      
00351     }
00352     delete [] dat;
00353   } 
00354 
00355   void CalculateForce(DataType& sum) {
00356     int me = MY_PROC; 
00357 
00358     double pi = 4.0*std::atan(1.0);
00359     
00360     // Read sphere info from F77 common block
00361     int Np;
00362     double xm, vm, cm, rd;
00363     f_combl_read(xm,vm,cm,rd,Np);
00364 
00365     // Calculate pressure
00366     for (register int l=0; l<=FineLevel(base::GH()); l++) {
00367       int Time = CurrentTime(base::GH(),l);
00368       int press_cnt = base::Dim()+4;      
00369       ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time, l, 
00370                                               press_cnt, base::t[l]);
00371     }
00372 
00373     int Npoints = Np*Np;
00374     DataType* p = new DataType[Npoints];
00375     point_type* xc = new point_type[Npoints];
00376 
00377     register int n;
00378 
00379     for (n=0; n<Np; n++) 
00380       for (register int m=0; m<Np; m++) {
00381         xc[n*Np+m](0) = xm; 
00382         xc[n*Np+m](1) = (2.*n/Np-1.)*rd;
00383         xc[n*Np+m](2) = (2.*m/Np-1.)*rd;
00384       }
00385     
00386     // Interpolate/extrapolate pressure values and assemble results on processor VizServer
00387     _Interpolation->PointsValues(base::Work(),base::t[0],Npoints,xc,p);
00388     _Interpolation->ArrayCombine(VizServer,Npoints,p);
00389 
00390     // Do the integration on the surface only on processor VizServer
00391     if (me == VizServer) {
00392       int vp = 0;
00393       sum = 0.;
00394       for (n=0; n<Npoints; n++)
00395         if (p[n]>-1.e37) {
00396           sum += p[n]; vp++;
00397         }
00398       sum /= vp;
00399       sum -= 101325.;
00400       sum *= pi*rd*rd;
00401 
00402       delete [] p;
00403       delete [] xc;
00404     }
00405 
00406 #ifdef DAGH_NO_MPI
00407 #else
00408     if (comm_service::dce() && comm_service::proc_num() > 1) {  
00409       int num = comm_service::proc_num();
00410       MPI_Status status;
00411       int R;
00412       int size = sizeof(DataType);
00413       int tag = 10001;
00414       if (me == VizServer) 
00415         for (int proc=0; proc<num; proc++) {
00416           if (proc==VizServer) continue;
00417           R = MPI_Send((void *)&sum, size, MPI_BYTE, proc, tag, comm_service::comm());
00418           if ( MPI_SUCCESS != R ) 
00419             comm_service::error_die( "SumCombine", "MPI_Send", R );
00420         }
00421       else {
00422         R = MPI_Recv((void *)&sum, size, MPI_BYTE, VizServer, tag, comm_service::comm(), &status);
00423         if ( MPI_SUCCESS != R ) 
00424           comm_service::error_die( "ArrayCombine", "MPI_Recv", R );
00425       }
00426     }
00427 #endif
00428   }
00429 
00430 protected:
00431   IntegratorSpecific _IntegratorSpecific;
00432   InitialConditionSpecific _InitialConditionSpecific;
00433   BoundaryConditionsSpecific _BoundaryConditionsSpecific;
00434   int OutputEvery, Nxc;
00435   point_type xc[OUTPUTPOINTSMAX];
00436   interpolation_type* _Interpolation; 
00437   ControlDevice SensorsCtrl;