vtf-logo

fsi/sfc-amroc/TubeCJBurnElastic/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@amroc.net
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 OWN_FLAGGING
00017 #define OWN_ELCGFMAMRSOLVER
00018 #include "ClpStdELCGFMProblem.h" 
00019 
00020 class FlaggingSpecific : 
00021   public AMRFlagging<VectorType,FixupType,FlagType,DIM> {
00022   typedef AMRFlagging<VectorType,FixupType,FlagType,DIM> base;
00023 public:
00024   FlaggingSpecific(base::solver_type& solver) : base(solver) {
00025       base::AddCriterion(new ByValue<VectorType,FlagType,DIM>());
00026       base::AddCriterion(new ScaledGradient<VectorType,FlagType,DIM>());
00027       base::AddCriterion(new LimiterType<VectorType,FlagType,DIM>());
00028       base::AddCriterion(new AbsoluteError<VectorType,FixupType,FlagType,DIM>(solver));
00029       base::AddCriterion(new RelativeError<VectorType,FixupType,FlagType,DIM>(solver));
00030 #ifdef f_flgout
00031       base::AddCriterion(new F77ByValue<VectorType,FlagType,DIM>(f_flgout));
00032       base::AddCriterion(new F77ScaledGradient<VectorType,FlagType,DIM>(f_flgout));
00033       base::AddCriterion(new F77LimiterType<VectorType,FlagType,DIM>(f_flgout));
00034       base::AddCriterion(new F77AbsoluteError<VectorType,FixupType,FlagType,DIM>(solver,f_flgout));
00035       base::AddCriterion(new F77RelativeError<VectorType,FixupType,FlagType,DIM>(solver,f_flgout));
00036 #endif
00037       for (int d=0; d<base::Dim(); d++) {
00038         dcmin[d] = 0; dcmax[d] = -1;
00039       }
00040     }
00041 
00042   ~FlaggingSpecific() { DeleteAllCriterions(); }
00043 
00044   virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00045     base::register_at(Ctrl, prefix);
00046     char VariableName[32];
00047     for (int d=0; d<base::Dim(); d++) {
00048       sprintf(VariableName,"DCellsMin(%d)",d+1);
00049       RegisterAt(base::LocCtrl,VariableName,dcmin[d]);
00050       sprintf(VariableName,"DCellsMax(%d)",d+1);
00051       RegisterAt(base::LocCtrl,VariableName,dcmax[d]);
00052     }
00053   }
00054   virtual void register_at(ControlDevice& Ctrl) {
00055     register_at(Ctrl, "");
00056   }
00057 
00058   virtual void SetFlags(const int Time, const int Level, double t, double dt) {
00059     base::SetFlags(Time, Level, t, dt);
00060     Coords lb(base::Dim(), dcmin), ub(base::Dim(), dcmax), s(base::Dim(),1);
00061     BBox bbox(lb*RefinedBy(base::GH(),MaxLevel(base::GH())),
00062               ub*RefinedBy(base::GH(),MaxLevel(base::GH())),
00063               s*StepSize(base::GH(),Level));
00064     forall (base::Flags(),Time,Level,c)  
00065       base::Flags()(Time,Level,c).equals(GoodPoint, bbox);
00066     end_forall    
00067   }
00068 
00069 private:
00070   int dcmin[DIM], dcmax[DIM];
00071 };
00072 
00073 
00074 template <class DataType, int dim>
00075 class F77CPTLevelSet : public CPTLevelSet<DataType,dim> {
00076   typedef CPTLevelSet<DataType,dim> base;
00077 public:
00078   typedef typename base::grid_fct_type grid_fct_type;  
00079   typedef typename base::grid_data_type grid_data_type;
00080   typedef generic_fortran_func generic_func_type; 
00081 
00082   typedef void (*lset_3_func_type) ( const INTEGER& maxmx, const INTEGER& maxmy, const INTEGER& maxmz,
00083                                      const INTEGER& mbc,
00084                                      const INTEGER& mx, const INTEGER& my, const INTEGER& mz,
00085                                      const DOUBLE x[], const DOUBLE y[], const DOUBLE z[],
00086                                      const DOUBLE& dx, const DOUBLE& dy, const DOUBLE& dz, 
00087                                      DataType q[], const DOUBLE& t);
00088 
00089   F77CPTLevelSet() : base(), f_lset(0) {}
00090   F77CPTLevelSet(generic_func_type lset) : base(), f_lset(lset) {}
00091 
00092   virtual void SetPhi(grid_fct_type& phi, const int Time, const int Level, double t) { 
00093     base::SetPhi(phi,Time,Level,t);
00094     forall (phi,Time,Level,c)
00095       SetGrid(phi(Time,Level,c),Level,t);    
00096     end_forall
00097   }
00098 
00099   virtual void SetGrid(grid_data_type& gdphi, const int& level, const double& t) {   
00100     assert(f_lset != (generic_func_type) 0);
00101     Coords ex = gdphi.extents();
00102     DCoords lbcorner = base::GH().worldCoords(gdphi.lower(), gdphi.stepsize());
00103     DCoords dx = base::GH().worldStep(gdphi.stepsize());
00104     int maxm[3], mx[3], d;
00105     DataType* x[3];
00106     for (d=0; d<dim; d++) {
00107       maxm[d] = ex(d);
00108       mx[d] = ex(d)-2*base::NGhosts();
00109       x[d] = new DataType[maxm[d]];
00110       for (int i=0; i<maxm[d]; i++) x[d][i] = (i+0.5)*dx(d)+lbcorner(d);
00111     }
00112 
00113     ((lset_3_func_type) f_lset)(AA(3,mx), base::NGhosts(), AA(3,mx), 
00114                                 AA(3,x), AA(3,dx), gdphi.data(), t);
00115 
00116     for (d=0; d<dim; d++) 
00117       delete [] x[d];
00118   }
00119 
00120   inline void SetFunc(generic_func_type lset) { f_lset = lset; }
00121   generic_func_type GetFunc() const { return f_lset; }
00122 
00123 protected:
00124   generic_func_type f_lset;
00125 };
00126 
00127 
00128 
00129 class FluidSolverSpecific : 
00130   public AMRELCGFMSolver<VectorType,FixupType,FlagType,DIM> {
00131   typedef VectorType::InternalDataType DataType;
00132   typedef AMRELCGFMSolver<VectorType,FixupType,FlagType,DIM> base;
00133 public:
00134   typedef F77GFMFileOutput<VectorType,FixupType,FlagType,DIM> output_type;
00135   typedef base::point_type point_type;
00136   typedef base::eul_comm_type eul_comm_type;
00137 
00138   FluidSolverSpecific() : base(_IntegratorSpecific, _InitialConditionSpecific, 
00139                                _BoundaryConditionsSpecific) {
00140     SetLevelTransfer(new F77LevelTransfer<VectorType,DIM>(f_prolong, f_restrict));
00141     SetFileOutput(new output_type(*this,f_flgout)); 
00142     SetFixup(new FixupSpecific(_IntegratorSpecific));
00143     SetFlagging(new FlaggingSpecific(*this)); 
00144     AddGFM(new GhostFluidMethod<VectorType,DIM>(
00145               new F77ELCGFMBoundary<VectorType,DIM>(f_ibndrfl,f_itrans,*this),
00146               new F77CPTLevelSet<DataType,DIM>(f_lset)));
00147     SetCoupleGFM(0);
00148     _ErasePressures = -1.e37;
00149     _MovePoints = -1.e37;
00150     _OffsetPoints = 0.;
00151   }  
00152  
00153   ~FluidSolverSpecific() {
00154     DeleteGFM(_GFM[0]);
00155     delete _Flagging;
00156     delete _Fixup;
00157   }
00158 
00159   virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00160     base::register_at(Ctrl,prefix);
00161     RegisterAt(base::LocCtrl,"ErasePressures",_ErasePressures);
00162     RegisterAt(base::LocCtrl,"MovePoints",_MovePoints);
00163     RegisterAt(base::LocCtrl,"OffsetPoints",_OffsetPoints);
00164   } 
00165   virtual void register_at(ControlDevice& Ctrl) {
00166     base::register_at(Ctrl);
00167   }
00168 
00169   virtual void SendBoundaryData() {
00170     START_WATCH
00171       for (register int l=0; l<=FineLevel(base::GH()); l++) 
00172         ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), 
00173                                                 CurrentTime(base::GH(),l), l, 
00174                                                 base::Dim()+4, base::t[l]);
00175     END_WATCH(FLUID_CPL_PRESSURE_CALCULATE)
00176     base::SendBoundaryData();
00177   }
00178 
00179   virtual void ModifyELCReceiveData(eul_comm_type* eulComm) {
00180     point_type *pos = (point_type *)(eulComm->getPositionsData());
00181     for (register int n=0; n<_eulComm->getNumberOfNodes(); n++)
00182       if (pos[n](2)<=_MovePoints)
00183         pos[n](2) += _OffsetPoints;   
00184   }
00185 
00186   virtual void ModifyELCSendData(eul_comm_type* eulComm) {
00187     double *press_data = eulComm->getPressuresData();
00188     const int* con = eulComm->getConnectivitiesData();
00189     const point_type *pos = reinterpret_cast<const point_type *>(eulComm->getPositionsData());
00190     register int n;
00191 
00192     if (base::_ThinStructure>0 || base::SendDataLocation()==0) {
00193       const point_type *centroid = 
00194         reinterpret_cast<const point_type *>(_eulComm->getFaceCentroidsData());
00195       for (n=0; n<eulComm->getNumberOfFaces(); n++)
00196         if (centroid[n](2)<=_ErasePressures)
00197           press_data[n] = std::numeric_limits<DataType>::max();
00198     }
00199     else 
00200       for (n=0; n<_eulComm->getNumberOfNodes(); n++)
00201         if (pos[n](2)<_ErasePressures)
00202           press_data[n] = std::numeric_limits<DataType>::max();
00203   }
00204 
00205 protected:
00206   IntegratorSpecific _IntegratorSpecific;
00207   InitialConditionSpecific _InitialConditionSpecific;
00208   BoundaryConditionsSpecific _BoundaryConditionsSpecific;