fsi/sfc-amroc/TubeCJBurnFlaps/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 "eulerznd3.h"
00010 
00011 #undef  NAUX
00012 #define NAUX  2
00013 
00014 #include "ClpProblem.h"
00015 
00016 #define MaxIntPoints (20)
00017 #define MaxDPoints (10)
00018 
00019 #define OWN_FLAGGING
00020 #define OWN_ELCGFMAMRSOLVER
00021 #include "ClpStdELCGFMProblem.h" 
00022 #include "AMRGFMInterpolation.h"
00023 
00024 class FlaggingSpecific : 
00025   public AMRFlagging<VectorType,FixupType,FlagType,DIM> {
00026   typedef AMRFlagging<VectorType,FixupType,FlagType,DIM> base;
00027 public:
00028   FlaggingSpecific(base::solver_type& solver) : base(solver) {
00029       base::AddCriterion(new ByValue<VectorType,FlagType,DIM>());
00030       base::AddCriterion(new ScaledGradient<VectorType,FlagType,DIM>());
00031       base::AddCriterion(new LimiterType<VectorType,FlagType,DIM>());
00032       base::AddCriterion(new AbsoluteError<VectorType,FixupType,FlagType,DIM>(solver));
00033       base::AddCriterion(new RelativeError<VectorType,FixupType,FlagType,DIM>(solver));
00034 #ifdef f_flgout
00035       base::AddCriterion(new F77ByValue<VectorType,FlagType,DIM>(f_flgout));
00036       base::AddCriterion(new F77ScaledGradient<VectorType,FlagType,DIM>(f_flgout));
00037       base::AddCriterion(new F77LimiterType<VectorType,FlagType,DIM>(f_flgout));
00038       base::AddCriterion(new F77AbsoluteError<VectorType,FixupType,FlagType,DIM>(solver,f_flgout));
00039       base::AddCriterion(new F77RelativeError<VectorType,FixupType,FlagType,DIM>(solver,f_flgout));
00040 #endif
00041       dcmin = new int[MaxDPoints*base::Dim()];
00042       dcmax = new int[MaxDPoints*base::Dim()];
00043       dclev = new int[MaxDPoints];
00044       for (int d=0; d<MaxDPoints*base::Dim(); d++) {
00045         dcmin[d] = 0; dcmax[d] = -1;
00046       }
00047       for (int i=0; i<MaxDPoints; i++)
00048         dclev[i] = 1000000;
00049     }
00050 
00051   ~FlaggingSpecific() { 
00052     DeleteAllCriterions(); 
00053     delete [] dcmin;
00054     delete [] dcmax;
00055   }
00056 
00057   virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00058     base::register_at(Ctrl, prefix);
00059     char VariableName[32];
00060     for (int i=0; i<MaxDPoints; i++) {
00061       for (int d=0; d<base::Dim(); d++) {
00062         sprintf(VariableName,"DCellsMin(%d,%d)",i+1,d+1);
00063         RegisterAt(base::LocCtrl,VariableName,dcmin[i*base::Dim()+d]);
00064         sprintf(VariableName,"DCellsMax(%d,%d)",i+1,d+1);
00065         RegisterAt(base::LocCtrl,VariableName,dcmax[i*base::Dim()+d]);
00066       }
00067       sprintf(VariableName,"DMinLevel(%d)",i+1);
00068       RegisterAt(base::LocCtrl,VariableName,dclev[i]);
00069     }
00070   }
00071   virtual void register_at(ControlDevice& Ctrl) {
00072     register_at(Ctrl, "");
00073   }
00074 
00075   virtual void SetFlags(const int Time, const int Level, double t, double dt) {
00076     base::SetFlags(Time, Level, t, dt);
00077     for (int i=0; i<MaxDPoints; i++) 
00078       if (Level>=dclev[i]) {
00079         Coords lb(base::Dim(), dcmin+i*base::Dim()), ub(base::Dim(), dcmax+i*base::Dim()), 
00080                   s(base::Dim(),1);
00081         BBox bbox(lb*RefinedBy(base::GH(),MaxLevel(base::GH())),
00082                   ub*RefinedBy(base::GH(),MaxLevel(base::GH())),
00083                   s*StepSize(base::GH(),Level));
00084         forall (base::Flags(),Time,Level,c)  
00085           base::Flags()(Time,Level,c).equals(GoodPoint, bbox);
00086         end_forall    
00087       }
00088   }
00089 
00090 private:
00091   int *dcmin, *dcmax, *dclev;
00092 };
00093 
00094 
00095 template <class DataType, int dim>
00096 class F77CPTLevelSet : public CPTLevelSet<DataType,dim> {
00097   typedef CPTLevelSet<DataType,dim> base;
00098 public:
00099   typedef typename base::grid_fct_type grid_fct_type;  
00100   typedef typename base::grid_data_type grid_data_type;
00101   typedef generic_fortran_func generic_func_type; 
00102 
00103   typedef void (*lset_3_func_type) ( const INTEGER& maxmx, const INTEGER& maxmy, const INTEGER& maxmz,
00104                                      const INTEGER& mbc,
00105                                      const INTEGER& mx, const INTEGER& my, const INTEGER& mz,
00106                                      const DOUBLE x[], const DOUBLE y[], const DOUBLE z[],
00107                                      const DOUBLE& dx, const DOUBLE& dy, const DOUBLE& dz, 
00108                                      DataType q[], const DOUBLE& t);
00109 
00110   F77CPTLevelSet() : base(), f_lset(0) {}
00111   F77CPTLevelSet(generic_func_type lset) : base(), f_lset(lset) {}
00112 
00113   virtual void SetPhi(grid_fct_type& phi, const int Time, const int Level, double t) { 
00114     base::SetPhi(phi,Time,Level,t);
00115     forall (phi,Time,Level,c)
00116       SetGrid(phi(Time,Level,c),Level,t);    
00117     end_forall
00118   }
00119 
00120   virtual void SetGrid(grid_data_type& gdphi, const int& level, const double& t) {   
00121     assert(f_lset != (generic_func_type) 0);
00122     Coords ex = gdphi.extents();
00123     DCoords lbcorner = base::GH().worldCoords(gdphi.lower(), gdphi.stepsize());
00124     DCoords dx = base::GH().worldStep(gdphi.stepsize());
00125     int maxm[3], mx[3], d;
00126     DataType* x[3];
00127     for (d=0; d<dim; d++) {
00128       maxm[d] = ex(d);
00129       mx[d] = ex(d)-2*base::NGhosts();
00130       x[d] = new DataType[maxm[d]];
00131       for (int i=0; i<maxm[d]; i++) x[d][i] = (i+0.5)*dx(d)+lbcorner(d);
00132     }
00133 
00134     ((lset_3_func_type) f_lset)(AA(3,mx), base::NGhosts(), AA(3,mx), 
00135                                 AA(3,x), AA(3,dx), gdphi.data(), t);
00136 
00137     for (d=0; d<dim; d++) 
00138       delete [] x[d];
00139   }
00140 
00141   inline void SetFunc(generic_func_type lset) { f_lset = lset; }
00142   generic_func_type GetFunc() const { return f_lset; }
00143 
00144 protected:
00145   generic_func_type f_lset;
00146 };
00147 
00148 
00149 
00150 class FluidSolverSpecific : 
00151   public AMRELCGFMSolver<VectorType,FixupType,FlagType,DIM> {
00152   typedef VectorType::InternalDataType DataType;
00153   typedef AMRELCGFMSolver<VectorType,FixupType,FlagType,DIM> base;
00154   typedef AMRGFMInterpolation<VectorType,FixupType,FlagType,DIM> interpolation_type;
00155 public:
00156   typedef F77GFMFileOutput<VectorType,FixupType,FlagType,DIM> output_type;
00157   typedef base::point_type point_type;
00158   typedef base::eul_comm_type eul_comm_type;
00159 
00160   FluidSolverSpecific() : base(_IntegratorSpecific, _InitialConditionSpecific, 
00161                                _BoundaryConditionsSpecific) {
00162     SetLevelTransfer(new F77LevelTransfer<VectorType,DIM>(f_prolong, f_restrict));
00163     SetFileOutput(new output_type(*this,f_flgout)); 
00164     SetFixup(new FixupSpecific(_IntegratorSpecific));
00165     SetFlagging(new FlaggingSpecific(*this)); 
00166     AddGFM(new GhostFluidMethod<VectorType,DIM>(
00167               new F77ELCGFMBoundary<VectorType,DIM>(f_ibndrfl,f_itrans,*this),
00168               new F77CPTLevelSet<DataType,DIM>(f_lset)));
00169     _Interpolation = new interpolation_type(*this);
00170     SetCoupleGFM(0);
00171     _ErasePressures = -1.e37;
00172     _MovePoints = -1.e37;
00173     _OffsetPoints = 0.;
00174     _NIntPoints = 0;
00175     _IntName = "ptrack.txt";
00176   }  
00177  
00178   ~FluidSolverSpecific() {
00179     DeleteGFM(_GFM[0]);
00180     delete _Flagging;
00181     delete _Fixup;
00182     delete _Interpolation;
00183   }
00184 
00185   virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00186     base::register_at(Ctrl,prefix);
00187     RegisterAt(base::LocCtrl,"ErasePressures",_ErasePressures);
00188     RegisterAt(base::LocCtrl,"MovePoints",_MovePoints);
00189     RegisterAt(base::LocCtrl,"OffsetPoints",_OffsetPoints);
00190     IntCtrl = base::LocCtrl.getSubDevice("TrackPressure");
00191     RegisterAt(IntCtrl,"NPoints",_NIntPoints);
00192     char VariableName[32];
00193     for (int nc=0; nc<MaxIntPoints; nc++) {
00194       for (int d=0; d<base::Dim(); d++) {
00195         std::sprintf(VariableName,"Point(%d,%d)",nc+1,d+1);
00196         RegisterAt(IntCtrl,VariableName,_IntPoints[nc](d)); 
00197         _IntPoints[nc](d) = 0.0;
00198       }
00199     }
00200     RegisterAt(IntCtrl,"FileName",_IntName);
00201   } 
00202   virtual void register_at(ControlDevice& Ctrl) {
00203     base::register_at(Ctrl);
00204   }
00205 
00206   virtual void SetupData() {
00207     base::SetupData();
00208     _Interpolation->SetupData(base::PGH(), base::NGhosts());
00209   }
00210 
00211   virtual void SendBoundaryData() {
00212     START_WATCH
00213       for (register int l=0; l<=FineLevel(base::GH()); l++) {
00214         int Time = CurrentTime(base::GH(),l);
00215         ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time,
00216                                                 l, base::Dim()+4, base::t[l]);
00217         if (CurrentTime(base::GH(),base::CouplingLevel) != Time) 
00218           ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time+TimeStep(base::U(),l),
00219                                                   l, base::Dim()+4, base::t[l]+base::dt[l]);      
00220       }
00221     END_WATCH(FLUID_CPL_PRESSURE_CALCULATE)
00222     base::SendBoundaryData();
00223   }
00224 
00225   virtual void ModifyELCReceiveData(eul_comm_type* eulComm) {
00226     point_type *pos = (point_type *)(eulComm->getPositionsData());
00227     for (register int n=0; n<_eulComm->getNumberOfNodes(); n++)
00228       if (pos[n](2)<=_MovePoints)
00229         pos[n](2) += _OffsetPoints;   
00230   }
00231 
00232   virtual void ModifyELCSendData(eul_comm_type* eulComm) {
00233     double *press_data = eulComm->getPressuresData();
00234     const int* con = eulComm->getConnectivitiesData();
00235     const point_type *pos = reinterpret_cast<const point_type *>(eulComm->getPositionsData());
00236     register int n;
00237 
00238     if (base::_ThinStructure>0 || base::SendDataLocation()==0) {
00239       const point_type *centroid = 
00240         reinterpret_cast<const point_type *>(_eulComm->getFaceCentroidsData());
00241       for (n=0; n<eulComm->getNumberOfFaces(); n++)
00242         if (centroid[n](2)<=_ErasePressures)
00243           press_data[n] = std::numeric_limits<DataType>::max();
00244     }
00245     else 
00246       for (n=0; n<_eulComm->getNumberOfNodes(); n++)
00247         if (pos[n](2)<_ErasePressures)
00248           press_data[n] = std::numeric_limits<DataType>::max();
00249   }
00250 
00251   virtual void Advance(double& t, double& dt) {
00252     base::Advance(t,dt);
00253     if (_NIntPoints<=0 || _NIntPoints>MaxIntPoints) return;
00254 
00255     int me = MY_PROC;
00256     // Calculate pressure
00257     for (register int l=0; l<=FineLevel(base::GH()); l++) {
00258       int Time = CurrentTime(base::GH(),l);
00259       int press_cnt = base::Dim()+4;
00260       ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time, l, 
00261                                               press_cnt, base::t[l]);
00262     }
00263 
00264     DataType* p = new DataType[_NIntPoints];
00265     _Interpolation->PointsValues(base::Work(),base::t[0],_NIntPoints,_IntPoints,p);
00266     _Interpolation->ArrayCombine(VizServer,_NIntPoints,p);
00267 
00268     // Append to output file only on processor VizServer
00269     if (me == VizServer) {
00270       std::ofstream outfile;
00271       std::ostream* out;
00272       outfile.open(_IntName.c_str(), std::ios::out | std::ios::app);
00273       out = new std::ostream(outfile.rdbuf());
00274       *out << base::t[0];
00275       for (int ns=0; ns<_NIntPoints; ns++)
00276         *out << " " << p[ns]; 
00277       *out << std::endl;
00278       outfile.close();
00279       delete out;      
00280     }
00281 
00282     delete [] p;
00283   } 
00284 
00285 protected:
00286   IntegratorSpecific _IntegratorSpecific;
00287   InitialConditionSpecific _InitialConditionSpecific;
00288   BoundaryConditionsSpecific _BoundaryConditionsSpecific;
00289   DataType _ErasePressures, _MovePoints, _OffsetPoints;
00290   interpolation_type* _Interpolation; 
00291   ControlDevice IntCtrl;
00292   int _NIntPoints;
00293   point_type _IntPoints[MaxIntPoints];
00294   std::string _IntName;