00001
00002
00003
00004
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 int Time = CurrentTime(base::GH(),l);
00173 ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time,
00174 l, base::Dim()+4, base::t[l]);
00175 if (CurrentTime(base::GH(),base::CouplingLevel) != Time)
00176 ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time+TimeStep(base::U(),l),
00177 l, base::Dim()+4, base::t[l]+base::dt[l]);
00178 }
00179 END_WATCH(FLUID_CPL_PRESSURE_CALCULATE)
00180 base::SendBoundaryData();
00181 }
00182
00183 virtual void ModifyELCReceiveData(eul_comm_type* eulComm) {
00184 point_type *pos = (point_type *)(eulComm->getPositionsData());
00185 for (register int n=0; n<_eulComm->getNumberOfNodes(); n++)
00186 if (pos[n](2)<=_MovePoints)
00187 pos[n](2) += _OffsetPoints;
00188 }
00189
00190 virtual void ModifyELCSendData(eul_comm_type* eulComm) {
00191 double *press_data = eulComm->getPressuresData();
00192 const int* con = eulComm->getConnectivitiesData();
00193 const point_type *pos = reinterpret_cast<const point_type *>(eulComm->getPositionsData());
00194 register int n;
00195
00196 if (base::_ThinStructure>0 || base::SendDataLocation()==0) {
00197 const point_type *centroid =
00198 reinterpret_cast<const point_type *>(_eulComm->getFaceCentroidsData());
00199 for (n=0; n<eulComm->getNumberOfFaces(); n++)
00200 if (centroid[n](2)<=_ErasePressures)
00201 press_data[n] = std::numeric_limits<DataType>::max();
00202 }
00203 else
00204 for (n=0; n<_eulComm->getNumberOfNodes(); n++)
00205 if (pos[n](2)<_ErasePressures)
00206 press_data[n] = std::numeric_limits<DataType>::max();
00207 }
00208
00209 protected:
00210 IntegratorSpecific _IntegratorSpecific;
00211 InitialConditionSpecific _InitialConditionSpecific;
00212 BoundaryConditionsSpecific _BoundaryConditionsSpecific;
00213 DataType _ErasePressures, _MovePoints, _OffsetPoints;