00001
00002
00003 #ifndef AMROC_FLUIDPROBLEM_H
00004 #define AMROC_FLUIDPROBLEM_H
00005 #define DIM 3
00006
00007 #include "LBMProblem.h"
00008 #include "LBMD3Q19.h"
00009
00010 typedef double DataType;
00011 typedef LBMD3Q19<DataType> LBMType;
00012
00013 #define OWN_LBMSCHEME
00014 #define OWN_AMRSOLVER
00015
00016
00017 #include "LBMStdProblem.h"
00018 #include "F77Interfaces/F77GFMBoundary.h"
00019 #include "AMRELCGFMSolver.h"
00020 #include "LBMGFMBoundary.h"
00021 #include "Interfaces/SchemeGFMFileOutput.h"
00022
00023 class LBMSpecific : public LBMType {
00024 typedef LBMType base;
00025 typedef ads::FixedArray<8,DataType> PVType;
00026 typedef ads::FixedArray<10,DataType> PVCType;
00027 public:
00028 LBMSpecific() : base(), omega(1.), u_p_avg(1.), l0_p(-1.) {
00029
00030 }
00031
00032 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00033 base::register_at(Ctrl,prefix);
00034 RegisterAt(base::LocCtrl,"Omega",omega);
00035 RegisterAt(base::LocCtrl,"u_p_avg",u_p_avg);
00036 RegisterAt(base::LocCtrl,"l0_p",l0_p);
00037
00038 }
00039
00040 virtual void SetupData(GridHierarchy* gh, const int& ghosts) {
00041 base::SetupData(gh,ghosts);
00042
00043 const double* bndry = base::GH().wholebndry();
00044 const int nbndry = base::GH().nbndry();
00045 if (l0_p <= 0. ) {
00046 l0_p=bndry[nbndry*(1+2*1)]-bndry[nbndry*(0+2*1)];
00047 if (nbndry>1) l0_p=bndry[nbndry*(1+2*1)]-bndry[1+nbndry*(1+2*1)];
00048
00049 }
00050 DataType omega = base::Omega(base::TimeScale());
00051
00052
00053 double Re_p=u_p_avg*l0_p/base::GasViscosity();
00054 double Re_lbm=(u_p_avg/base::VelocityScale())*(l0_p/base::LengthScale())/base::LatticeViscosity(omega);
00055 ( comm_service::log() << "l0_p=" << l0_p << " u_p_avg=" << u_p_avg
00056 << " Re_p=" << Re_p << " Re_lbm=" << Re_lbm
00057 << "\n omega=" << omega
00058 << " SpeedUp=" << base::SpeedUp()
00059 << " u_avg_LBM=" << u_p_avg/base::VelocityScale()
00060 << " SpeedRatio=" << u_p_avg/base::VelocityScale()*std::sqrt(DataType(3.))
00061 << " base::TO()=" << base::T0()
00062 << "\n Re_p-Re_lbm=" << Re_p-Re_lbm
00063 << std::endl ).flush();
00064
00065 if ( std::fabs(Re_p-Re_lbm)>DBL_EPSILON*LB_FACTOR*LB_FACTOR*LB_FACTOR ) {
00066 std::cerr << "Physical scaling requires fabs(Re_p-Re_lbm)<DBL_EPSILON*LB_FACTOR*LB_FACTOR*LB_FACTOR !" << std::endl;
00067 std::cerr << "Physical scaling requires fabs(Re_p-Re_lbm)<DBL_EPSILON*LB_FACTOR*LB_FACTOR*LB_FACTOR !" << std::endl;
00068 std::cerr << "Re_p " << Re_p << " Re_lbm " << Re_lbm << " std::fabs(Re_p-Re_lbm) " << std::fabs(Re_p-Re_lbm) << std::endl;
00069 std::exit(-1);
00070 }
00071 if (omega<0.5 || omega>2.) {
00072 std::cerr << "0.5 <= omega <= 2 required." << std::endl;
00073 std::exit(-1);
00074 }
00075 base::SetGas(base::GasDensity(),base::GasViscosity(),base::GasSpeedofSound());
00076 ( comm_service::log() << "FluidProblem SetupData() D3Q19: Gas_rho=" << base::GasDensity() << " Gas_Cs=" << base::GasSpeedofSound()
00077 << " Gas_nu=" << base::GasViscosity() << std::endl ).flush();
00078
00079 }
00080
00081 virtual double Step(vec_grid_data_type &fvec, vec_grid_data_type &ovec, vec_grid_data_type* Flux[],
00082 const double& t, const double& dt, const int& mpass) const {
00083
00084 LBMType::Step(fvec, ovec, Flux, t, dt, mpass);
00085 return DataType(1.);
00086 }
00087
00088 private:
00089 DataType omega, u_p_avg, l0_p;
00090
00091 protected:
00092
00093 };
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286 template <class LBMType, int dim>
00287 class LBMELCGFMBoundary : public LBMGFMBoundary<LBMType,dim> {
00288 typedef typename LBMType::MicroType MicroType;
00289 typedef typename MicroType::InternalDataType DataType;
00290 typedef LBMGFMBoundary<LBMType,dim> base;
00291 public:
00292 typedef typename base::vec_grid_data_type vec_grid_data_type;
00293 typedef typename base::grid_data_type grid_data_type;
00294 typedef typename base::point_type point_type;
00295 typedef AMRELCGFMSolver<MicroType,FixupType,FlagType,dim> amr_elc_solver;
00296
00297 LBMELCGFMBoundary(LBMType &scheme, amr_elc_solver& solver) : base(scheme), _solver(solver) {
00298 base::SetNAux(base::Dim());
00299 }
00300
00301 virtual void SetBndryAux(vec_grid_data_type& gdu, grid_data_type& gdphi,
00302 const MicroType* u, DataType* aux, const int& Level,
00303 double t, const int& nc, const int* idx,
00304 const point_type* xc, const DataType* distance,
00305 const point_type* normal) {
00306 _solver.ExtractBoundaryVelocities(gdu.bbox(),nc,xc,reinterpret_cast<point_type *>(aux));
00307 }
00308
00309 protected:
00310 amr_elc_solver& _solver;
00311 };
00312
00313 class SolverObjects {
00314 public:
00315 SolverObjects() {
00316 _LBMSpecific = new LBMSpecific();
00317 _IntegratorSpecific = new IntegratorSpecific(*_LBMSpecific);
00318 _InitialConditionSpecific = new InitialConditionSpecific(*_LBMSpecific);
00319 _BoundaryConditionsSpecific = new BoundaryConditionsSpecific(*_LBMSpecific);
00320 }
00321
00322 ~SolverObjects() {
00323 delete _BoundaryConditionsSpecific;
00324 delete _InitialConditionSpecific;
00325 delete _IntegratorSpecific;
00326 delete _LBMSpecific;
00327 }
00328
00329 protected:
00330 LBMSpecific *_LBMSpecific;
00331 IntegratorSpecific *_IntegratorSpecific;
00332 InitialConditionSpecific *_InitialConditionSpecific;
00333 BoundaryConditionsSpecific *_BoundaryConditionsSpecific;
00334 };
00335
00336 class FluidSolverSpecific : protected SolverObjects,
00337 public AMRELCGFMSolver<MicroType,FixupType,FlagType,DIM> {
00338 typedef MicroType::InternalDataType DataType;
00339 typedef AMRELCGFMSolver<MicroType,FixupType,FlagType,DIM> base;
00340 typedef AMRInterpolation<MicroType,DIM> interpolation_type;
00341 typedef interpolation_type::point_type point_type;
00342 typedef SchemeGFMFileOutput<LBMType,FixupType,FlagType,DIM> output_type;
00343 typedef ads::FixedArray<10,DataType> PVCType;
00344 public:
00345 FluidSolverSpecific() : SolverObjects(), base(*SolverObjects::_IntegratorSpecific,
00346 *SolverObjects::_InitialConditionSpecific,
00347 *SolverObjects::_BoundaryConditionsSpecific) {
00348 SetLevelTransfer(new LBMF77LevelTransfer<LBMType,DIM>(SolverObjects::_IntegratorSpecific->Scheme(), f_prolong, f_restrict));
00349 SetFileOutput(new output_type(*this,SolverObjects::_IntegratorSpecific->Scheme()));
00350
00351 SetFixup(new FixupSpecific(SolverObjects::_IntegratorSpecific->Scheme()));
00352 SetFlagging(new FlaggingSpecific(*this,SolverObjects::_IntegratorSpecific->Scheme()));
00353 AddGFM(new GhostFluidMethod<MicroType,DIM>(new LBMELCGFMBoundary<LBMType,DIM>(SolverObjects::_IntegratorSpecific->LBM(),*this),
00354 new CPTLevelSet<DataType,DIM>()));
00355
00356 SetCoupleGFM(0);
00357
00358 }
00359
00360 ~FluidSolverSpecific() {
00361 DeleteGFM(_GFM[0]);
00362 delete _LevelTransfer;
00363 delete _Flagging;
00364 delete _Fixup;
00365 delete _FileOutput;
00366 delete _BoundaryConditionsSpecific;
00367 _BoundaryConditionsSpecific = (BoundaryConditionsSpecific *) 0;
00368 delete _InitialConditionSpecific;
00369 _InitialConditionSpecific = (InitialConditionSpecific *) 0;
00370 delete _IntegratorSpecific;
00371 _IntegratorSpecific = (IntegratorSpecific *) 0;
00372 delete _Interpolation;
00373 _Interpolation = (AMRELCGFMSolver<MicroType,FixupType,FlagType,DIM>::interpolation_type *) 0;
00374 }
00375
00376 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00377 base::register_at(Ctrl,prefix);
00378
00379 }
00380 virtual void register_at(ControlDevice& Ctrl) {
00381 base::register_at(Ctrl);
00382 }
00383
00384 virtual void SendBoundaryData() {
00385 START_WATCH
00386 for (register int l=0; l<=FineLevel(base::GH()); l++) {
00387 int Time = CurrentTime(base::GH(),l);
00388 ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time,
00389 l, base::Dim()+4, base::t[l]);
00390 if (CurrentTime(base::GH(),base::CouplingLevel) != Time)
00391 ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time+TimeStep(base::U(),l),
00392 l, base::Dim()+4, base::t[l]+base::dt[l]);
00393 }
00394 END_WATCH(FLUID_CPL_PRESSURE_CALCULATE)
00395 base::SendBoundaryData();
00396 }
00397
00398 virtual void SetupData() {
00399 ( comm_service::log() << " fluid setup begun" << std::endl ).flush();
00400 base::SetupData();
00401 base::AdaptBndTimeInterpolate = 0;
00402 base::NAMRTimeSteps = 1;
00403 base::Step[0].LastTime = 1.e37;
00404 base::Step[0].VariableTimeStepping = -1;
00405 base::Step[0].dtv[0] = ((LBMIntegrator<LBMType,DIM> &)Integrator_()).LBM().TimeScale();
00406
00407 _Interpolation->SetupData(base::PGH(), base::NGhosts());
00408
00409 }
00410
00411 protected:
00412
00413 };