fsi/sfc-amroc/WaterBlastFracture/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 CPTLevelSet<DataType,DIM>()));
00178     AddGFM(new GhostFluidMethod<VectorType,DIM>(
00179               new F77GFMBoundary<VectorType,DIM>(f_ibndrfl,f_itrans),
00180               new F77GFMLevelSet<DataType,DIM>(f_lset)));
00181     AddGFM(new GhostFluidMethod<VectorType,DIM>(
00182               new F77GFMBoundary<VectorType,DIM>(f_ibndrfl,f_itrans,f_ipvel_piston),
00183               new F77GFMLevelSet<DataType,DIM>(f_lset_piston)));
00184     SetCoupleGFM(0);
00185     _Interpolation = new interpolation_type(*this);
00186   }  
00187  
00188   ~FluidSolverSpecific() {
00189     DeleteGFM(_GFM[2]);
00190     DeleteGFM(_GFM[1]);
00191     DeleteGFM(_GFM[0]);
00192     delete _Flagging;
00193     delete _Fixup;
00194     delete _Interpolation;
00195   }
00196 
00197   virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00198     base::register_at(Ctrl,prefix);
00199     SensorsCtrl = base::LocCtrl.getSubDevice("SensorsOutput");
00200     RegisterAt(SensorsCtrl,"OutputEvery",OutputEvery);
00201     RegisterAt(SensorsCtrl,"NSensors",Nxc);
00202     char VariableName[32];
00203     for (register int i=0; i<OUTPUTPOINTSMAX; i++) 
00204       for (register int d=0; d<base::Dim(); d++) {
00205         std::sprintf(VariableName,"xc(%d,%d)",i+1,d+1);
00206         RegisterAt(SensorsCtrl,VariableName,xc[i](d));
00207       }
00208   } 
00209   virtual void register_at(ControlDevice& Ctrl) {
00210     base::register_at(Ctrl);
00211   }
00212 
00213   virtual void SetupData() {
00214     base::SetupData();
00215     _Interpolation->SetupData(base::PGH(), base::NGhosts());
00216     if (OutputEvery<0) OutputEvery=0;
00217     if (Nxc<0) Nxc=0;
00218     if (Nxc>OUTPUTPOINTSMAX) Nxc=OUTPUTPOINTSMAX;
00219   }
00220 
00221   virtual void Restart(double& t, double& dt) {
00222     std::ifstream ifs("boundary.cp", std::ios::in);
00223     int Np;
00224     double xm, vm, cm, rd;
00225     f_combl_read(xm,vm,cm,rd,Np);
00226     ifs >> xm >> vm;
00227     ifs.close();
00228     f_combl_write(xm,vm,cm,rd,Np);
00229     base::Restart(t,dt);
00230   }
00231   virtual void Checkpointing() {
00232     int me = MY_PROC; 
00233     if (me == VizServer) {
00234       int Np;
00235       double xm, vm, cm, rd;
00236       f_combl_read(xm,vm,cm,rd,Np);
00237       std::ofstream ofs("boundary.cp", std::ios::out);
00238       ofs << xm << " " << vm << std::endl; 
00239       ofs.close();
00240     }
00241     base::Checkpointing();
00242   }
00243 
00244   virtual void SendBoundaryData() {
00245     START_WATCH
00246       for (register int l=0; l<=FineLevel(base::GH()); l++) {
00247         int Time = CurrentTime(base::GH(),l);
00248         ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time,
00249                                                 l, base::Dim()+4, base::t[l]);
00250         if (CurrentTime(base::GH(),base::CouplingLevel) != Time) 
00251           ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time+TimeStep(base::U(),l),
00252                                                   l, base::Dim()+4, base::t[l]+base::dt[l]);      
00253       }
00254     END_WATCH(FLUID_CPL_PRESSURE_CALCULATE)
00255     base::SendBoundaryData();
00256   }
00257 
00258   virtual void ModifyELCSendData(eul_comm_type* eulComm) {
00259     DataType distance = -GFM(CoupleGFM()).Boundary().Cutoff();
00260     int Npoints = _eulComm->getNumberOfFaces();
00261     double *press_data = _eulComm->getPressuresData();
00262     const int* con = _eulComm->getConnectivitiesData();
00263     const point_type *pos = reinterpret_cast<const point_type *>(_eulComm->getPositionsData());
00264     for (register int n=0; n<Npoints; n++) {
00265       bool inside = true;
00266       for (register int d=0; d<base::Dim(); d++) {
00267         const point_type &p = pos[con[base::Dim()*n+d]]; 
00268         if (std::sqrt(p(1)*p(1)+p(2)*p(2))>0.0323) {
00269           inside = false;
00270           break;
00271         }
00272       }
00273       if (!inside) 
00274         press_data[n] = std::numeric_limits<DataType>::max();
00275     }
00276   }
00277   
00278   // Overloading of Advance() to add pointwise output after each level 0 time step
00279   // and propagate piston
00280   virtual void Advance(double& t, double& dt) {
00281 
00282     int me = MY_PROC; 
00283     DataType force_old;
00284     if (base::GFM(2).IsUsed())
00285       CalculateForce(force_old);
00286 
00287     base::Advance(t,dt);
00288 
00289     if (base::GFM(2).IsUsed()) {
00290       DataType a, force_new;
00291       CalculateForce(force_new);
00292       
00293       int Np;
00294       double xm, vm, cm, rd;
00295       f_combl_read(xm,vm,cm,rd,Np);
00296       
00297       a = 0.5/cm*(force_old+force_new);
00298       vm += a*dt;
00299       xm += vm*dt;
00300       
00301       f_combl_write(xm,vm,cm,rd,Np);
00302     
00303       // Append to output file only on processor VizServer
00304       if (me == VizServer) {
00305         std::ofstream outfile;
00306         std::ostream* out;
00307         std::string fname = "boundary.txt";
00308         outfile.open(fname.c_str(), std::ios::out | std::ios::app);
00309         out = new std::ostream(outfile.rdbuf());
00310         *out << t << " " << force_new << " " << xm << " " << vm << std::endl;
00311         outfile.close();
00312         delete out;      
00313       }
00314     }
00315 
00316     if (!OutputEvery) return;
00317 
00318     int TimeZero  = CurrentTime(base::GH(),0);
00319     TimeZero /= RefinedBy(base::GH(),MaxLevel(base::GH()));
00320     if (TimeZero % OutputEvery != 0) return;
00321 
00322     VectorType* dat = new VectorType[Nxc];
00323 
00324     // Interpolate/extrapolate pressure values and assemble results on processor VizServer
00325     _Interpolation->PointsValues(base::U(),t,Nxc,xc,dat);
00326     int NValues = Nxc*base::NEquations();
00327     _Interpolation->ArrayCombine(VizServer,NValues,&(dat[0](0)));
00328 
00329     // Append to output file only on processor VizServer
00330     if (MY_PROC == VizServer) {
00331       DataType* output_data = new DataType[_FileOutput->Ncnt()*Nxc];
00332       BBox bb(3,0,0,0,Nxc-1,0,0,1,1,1);
00333       vec_grid_data_type val(bb,dat);
00334       for (int cnt=1;cnt<=_FileOutput->Ncnt(); cnt++) {
00335         if (_FileOutput->OutputName(cnt-1).c_str()[0] == '-') continue;
00336         grid_data_type output(bb,output_data+(cnt-1)*Nxc);
00337         ((output_type::out_3_func_type) (((output_type*) base::_FileOutput)->GetOutFunc()))
00338           (FA(3,val),FA(3,output),BOUNDING_BOX(bb),base::NEquations(),cnt,t);
00339         DataType* dummy;
00340         output.deallocate(dummy);
00341       }
00342       VectorType* dummy;
00343       val.deallocate(dummy);
00344       std::ofstream outfile;
00345       std::ostream* out;
00346       std::string fname = "sensors.txt";
00347       outfile.open(fname.c_str(), std::ios::out | std::ios::app);
00348       out = new std::ostream(outfile.rdbuf());
00349       *out << t << " ";
00350       for (register int j=0;j<Nxc; j++)
00351         for (int cnt=0;cnt<_FileOutput->Ncnt(); cnt++) 
00352           if (_FileOutput->OutputName(cnt).c_str()[0] != '-') 
00353             *out << " " << output_data[cnt*Nxc+j];
00354       *out << std::endl;
00355       delete [] output_data;
00356       outfile.close();
00357       delete out;      
00358     }
00359     delete [] dat;
00360   } 
00361 
00362   void CalculateForce(DataType& sum) {
00363     int me = MY_PROC; 
00364 
00365     double pi = 4.0*std::atan(1.0);
00366     
00367     // Read sphere info from F77 common block
00368     int Np;
00369     double xm, vm, cm, rd;
00370     f_combl_read(xm,vm,cm,rd,Np);
00371 
00372     // Calculate pressure
00373     for (register int l=0; l<=FineLevel(base::GH()); l++) {
00374       int Time = CurrentTime(base::GH(),l);
00375       int press_cnt = base::Dim()+4;      
00376       ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time, l, 
00377                                               press_cnt, base::t[l]);
00378     }
00379 
00380     int Npoints = Np*Np;
00381     DataType* p = new DataType[Npoints];
00382     point_type* xc = new point_type[Npoints];
00383 
00384     register int n;
00385 
00386     for (n=0; n<Np; n++) 
00387       for (register int m=0; m<Np; m++) {
00388         xc[n*Np+m](0) = xm; 
00389         xc[n*Np+m](1) = (2.*n/Np-1.)*rd;
00390         xc[n*Np+m](2) = (2.*m/Np-1.)*rd;
00391       }
00392     
00393     // Interpolate/extrapolate pressure values and assemble results on processor VizServer
00394     _Interpolation->PointsValues(base::Work(),base::t[0],Npoints,xc,p);
00395     _Interpolation->ArrayCombine(VizServer,Npoints,p);
00396 
00397     // Do the integration on the surface only on processor VizServer
00398     if (me == VizServer) {
00399       int vp = 0;
00400       sum = 0.;
00401       for (n=0; n<Npoints; n++)
00402         if (p[n]>-1.e37) {
00403           sum += p[n]; vp++;
00404         }
00405       sum /= vp;
00406       sum -= 101325.;
00407       sum *= pi*rd*rd;
00408 
00409       delete [] p;
00410       delete [] xc;
00411     }
00412 
00413 #ifdef DAGH_NO_MPI
00414 #else
00415     if (comm_service::dce() && comm_service::proc_num() > 1) {  
00416       int num = comm_service::proc_num();
00417       MPI_Status status;
00418       int R;
00419       int size = sizeof(DataType);
00420       int tag = 10001;
00421       if (me == VizServer) 
00422         for (int proc=0; proc<num; proc++) {
00423           if (proc==VizServer) continue;
00424           R = MPI_Send((void *)&sum, size, MPI_BYTE, proc, tag, comm_service::comm());
00425           if ( MPI_SUCCESS != R ) 
00426             comm_service::error_die( "SumCombine", "MPI_Send", R );
00427         }
00428       else {
00429         R = MPI_Recv((void *)&sum, size, MPI_BYTE, VizServer, tag, comm_service::comm(), &status);
00430         if ( MPI_SUCCESS != R ) 
00431           comm_service::error_die( "ArrayCombine", "MPI_Recv", R );
00432       }
00433     }
00434 #endif
00435   }
00436 
00437 protected:
00438   IntegratorSpecific _IntegratorSpecific;
00439   InitialConditionSpecific _InitialConditionSpecific;
00440   BoundaryConditionsSpecific _BoundaryConditionsSpecific;
00441   int OutputEvery, Nxc;
00442   point_type xc[OUTPUTPOINTSMAX];
00443   interpolation_type* _Interpolation; 
00444   ControlDevice SensorsCtrl;