00001
00002
00003 #ifndef AMROC_FLUIDPROBLEM_H
00004 #define AMROC_FLUIDPROBLEM_H
00005 #define DIM 3
00006
00007 #ifdef DEBUG_PRINT_ELC_VALUES
00008 #undef DEBUG_PRINT_ELC_VALUES
00009 #endif
00010
00011 #undef DEBUG_PRINT_ELC_VALUES
00012
00013
00014 #define NEQUATIONS 8 // LBM equations for gases in 3D (5 fields),
00015 #define NVARS 5 // Total number of active variables
00016 #define SJSON
00017
00018 #include "LBMProblem.h"
00019 #include "LBMD3Q19.h"
00020
00021 typedef double DataType;
00022 typedef LBMD3Q19<DataType> LBMType;
00023
00024 #define OWN_LBMSCHEME
00025 #define OWN_AMRSOLVER
00026 #define OWN_INITIALCONDITION
00027 #define OWN_BOUNDARYCONDITION
00028 #include "LBMStdProblem.h"
00029 #include "F77Interfaces/F77GFMBoundary.h"
00030 #include "AMRELCGFMSolver.h"
00031 #include "LBMGFMBoundary.h"
00032 #include "Interfaces/SchemeGFMFileOutput.h"
00033
00034 #include "AMRInterpolation.h"
00035 #include "AMRStatistics.h"
00036
00037 class LBMSpecific : public LBMType {
00038 typedef LBMType base;
00039 typedef ads::FixedArray<8,DataType> PVType;
00040 typedef ads::FixedArray<10,DataType> PVCType;
00041 typedef MicroType::InternalDataType DataType;
00042 typedef AMRInterpolation<MicroType,DIM> interpolation_type;
00043 public:
00044 LBMSpecific() : base(), omega(1.), u_p_avg(1.), l0_p(-1.) {}
00045
00046 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00047 base::register_at(Ctrl,prefix);
00048 RegisterAt(base::LocCtrl,"Omega",omega);
00049 RegisterAt(base::LocCtrl,"u_p_avg",u_p_avg);
00050 RegisterAt(base::LocCtrl,"l0_p",l0_p);
00051
00052 }
00053
00054 virtual void SetupData(GridHierarchy* gh, const int& ghosts) {
00055 base::SetupData(gh,ghosts);
00056 const double* bndry = base::GH().wholebndry();
00057 const int nbndry = base::GH().nbndry();
00058 if (l0_p <= 0. ) {
00059 l0_p=bndry[nbndry*(1+2*1)]-bndry[nbndry*(0+2*1)];
00060 if (nbndry>1) l0_p=bndry[nbndry*(1+2*1)]-bndry[1+nbndry*(1+2*1)];
00061 }
00062 DataType omega = base::Omega(base::TimeScale());
00063 double Re_p=u_p_avg*l0_p/base::GasViscosity();
00064 double Re_lbm=(u_p_avg/base::VelocityScale())*(l0_p/base::LengthScale())/base::LatticeViscosity(omega);
00065 std::cout << "l0_p=" << l0_p << " u_p_avg=" << u_p_avg
00066 << " Re_p=" << Re_p << " Re_lbm=" << Re_lbm
00067 << "\n omega=" << omega
00068 << " SpeedUp=" << base::SpeedUp()
00069 << " u_avg_LBM=" << u_p_avg/base::VelocityScale()
00070 << " SpeedRatio=" << u_p_avg/base::VelocityScale()*std::sqrt(DataType(3.))
00071 << " base::TO()=" << base::T0()
00072 << "\n Re_p-Re_lbm=" << Re_p-Re_lbm
00073 << std::endl;
00074 if ( std::fabs(Re_p-Re_lbm)>DBL_EPSILON*LB_FACTOR*LB_FACTOR*LB_FACTOR ) {
00075 std::cerr << "Physical scaling requires fabs(Re_p-Re_lbm)<DBL_EPSILON*LB_FACTOR*LB_FACTOR*LB_FACTOR !" << std::endl;
00076 std::cerr << "Physical scaling requires fabs(Re_p-Re_lbm)<DBL_EPSILON*LB_FACTOR*LB_FACTOR*LB_FACTOR !" << std::endl;
00077 std::cerr << "Re_p " << Re_p << " Re_lbm " << Re_lbm << " std::fabs(Re_p-Re_lbm) " << std::fabs(Re_p-Re_lbm) << std::endl;
00078 std::exit(-1);
00079 }
00080 if (omega<0.5 || omega>2.) {
00081 std::cerr << "0.5 <= omega <= 2 required." << std::endl;
00082 std::exit(-1);
00083 }
00084 base::SetGas(base::GasDensity(),base::GasViscosity(),base::GasSpeedofSound());
00085 std::cout << "FluidProblem SetupData() D3Q19: Gas_rho=" << base::GasDensity() << " Gas_Cs=" << base::GasSpeedofSound()
00086 << " Gas_nu=" << base::GasViscosity() << std::endl;
00087
00088 }
00089
00090 virtual double Step(vec_grid_data_type &fvec, vec_grid_data_type &ovec, vec_grid_data_type* Flux[],
00091 const double& t, const double& dt, const int& mpass) const {
00092
00093 LBMType::Step(fvec, ovec, Flux, t, dt, mpass);
00094 return DataType(1.);
00095 }
00096
00097
00098
00099 private:
00100 DataType omega, u_p_avg, l0_p;
00101
00102 protected:
00103 int track_every;
00104 };
00105
00106 class InitialConditionSpecific : public LBMInitialCondition<LBMType,DIM> {
00107 typedef LBMInitialCondition<LBMType,DIM> base;
00108 public:
00109 InitialConditionSpecific(LBMType &lbm) : base(lbm) {}
00110 typedef base::MicroType MicroType;
00111 typedef base::MacroType MacroType;
00112
00113 virtual void SetGrid(vec_grid_data_type& fvec, grid_data_type& workvec, const int& level) {
00114 MacroType q, q_;
00115 DataType elevation = 15., profile, height, ldx;
00116 BBox wb = base::GH().wholebbox();
00117 BBox bb = fvec.bbox();
00118 DCoords wc;
00119
00120 q_(0) = aux[0]/LBM().DensityScale();
00121 q_(1) = aux[1]/LBM().VelocityScale();
00122 q_(2) = aux[2]/LBM().VelocityScale();
00123 q_(3) = aux[3]/LBM().VelocityScale();
00124
00125 MicroType *f = (MicroType *)fvec.databuffer();
00126
00127 int Nx = fvec.extents(0), Ny = fvec.extents(1);
00128 int b = fvec.bottom();
00129 assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1));
00130 int mxs = std::max(bb.lower(0)/bb.stepsize(0),fvec.lower(0)/fvec.stepsize(0));
00131 int mys = std::max(bb.lower(1)/bb.stepsize(1),fvec.lower(1)/fvec.stepsize(1));
00132 int mzs = std::max(bb.lower(2)/bb.stepsize(2),fvec.lower(2)/fvec.stepsize(2));
00133 int mxe = std::min(bb.upper(0)/bb.stepsize(0),fvec.upper(0)/fvec.stepsize(0));
00134 int mye = std::min(bb.upper(1)/bb.stepsize(1),fvec.upper(1)/fvec.stepsize(1));
00135 int mze = std::min(bb.upper(2)/bb.stepsize(2),fvec.upper(2)/fvec.stepsize(2));
00136 int lc_indx[3];
00137
00138 lc_indx[0] = base::GH().worldCoords(fvec.lower(),fvec.stepsize())(0);
00139 lc_indx[1] = base::GH().worldCoords(fvec.lower(),fvec.stepsize())(1);
00140
00141 ldx = DataType(fvec.stepsize(2))/DataType(wb.stepsize(2))*LBM().L0();
00142
00143 for (register int k=mzs; k<=mze; k++) {
00144 lc_indx[2] = (int) std::floor(DataType(k)*DataType(ldx));
00145 wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
00146 height = DataType(wc(2))*DataType(wb.stepsize(2));
00147
00148
00149 if ( height < elevation ) {
00150 profile = DataType(height*height)/DataType(elevation*elevation);
00151
00152 q(0) = q_(0);
00153 q(1) = q_(1)*profile;
00154 q(2) = q_(2)*profile;
00155 q(3) = q_(3)*profile;
00156 }
00157 else {
00158 q = q_;
00159
00160 }
00161 for (register int j=mys; j<=mye; j++)
00162 for (register int i=mxs; i<=mxe; i++)
00163 f[b+LBM().idx(i,j,k,Nx,Ny)] = LBM().Equilibrium(q);
00164 }
00165
00166 }
00167 };
00168
00169 class BoundaryConditionsSpecific : public LBMBoundaryConditions<LBMType,DIM> {
00170 typedef LBMBoundaryConditions<LBMType,DIM> base;
00171 public:
00172 typedef base::MicroType MicroType;
00173 typedef base::MacroType MacroType;
00174 BoundaryConditionsSpecific(LBMType &lbm) : base(lbm) {}
00175
00176 DataType caux[6];
00177 int ctype, cside, cscaling, ncaux;
00178
00179 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00180 base::register_at(Ctrl, prefix);
00181 base::LocCtrl = Ctrl.getSubDevice(prefix+"BoundaryConditionCustom");
00182 RegisterAt(base::LocCtrl,"Type",ctype);
00183 RegisterAt(base::LocCtrl,"Side",cside);
00184 RegisterAt(base::LocCtrl,"Scaling",cscaling);
00185 RegisterAt(base::LocCtrl,"NCAux",ncaux);
00186 char VariableName[16];
00187 for (int i=0; i<6; i++) {
00188 caux[i]=DataType(1.);
00189 std::sprintf(VariableName,"CAux(%d)",i+1);
00190 RegisterAt(base::LocCtrl,VariableName,caux[i]);
00191 }
00192 }
00193
00194 virtual void SetupData(GridHierarchy* gh, const int& ghosts) {
00195 base::SetupData(gh,ghosts);
00196 const double* bndry = base::GH().wholebndry();
00197 const int nbndry = base::GH().nbndry();
00198 DataType l0_p=bndry[nbndry*(1+2*1)]-bndry[nbndry*(0+2*1)];
00199 if (nbndry>1) l0_p=bndry[nbndry*(1+2*1)]-bndry[1+nbndry*(1+2*1)];
00200 DataType omega = LBM().Omega(LBM().TimeScale());
00201
00202 DataType u_p_avg=caux[1];
00203 double Re_p=u_p_avg*l0_p/LBM().GasViscosity();
00204 double Re_lbm=(u_p_avg/LBM().VelocityScale())*(l0_p/LBM().LengthScale())/LBM().LatticeViscosity(omega);
00205 ( comm_service::log() << "l0_p=" << l0_p << " u_p_avg=" << u_p_avg
00206 << " Re_p=" << Re_p << " Re_lbm=" << Re_lbm
00207 << "\n omega=" << omega
00208 << " SpeedUp=" << LBM().SpeedUp()
00209 << " u_avg_LBM=" << u_p_avg/LBM().VelocityScale()
00210 << " SpeedRatio=" << u_p_avg/LBM().VelocityScale()*std::sqrt(DataType(3.))
00211 << " LBM().TO()=" << LBM().T0()
00212 << "\n Re_p-Re_lbm=" << Re_p-Re_lbm
00213 << std::endl ).flush();
00214 assert(std::fabs(Re_p-Re_lbm)<DBL_EPSILON*LB_FACTOR*LB_FACTOR*LB_FACTOR);
00215 }
00216
00217 virtual void SetBndry(vec_grid_data_type &fvec, const int& level, const BBox &bb,
00218 const int &dir, const double& time) {
00219
00220
00221
00222 base::SetBndry(fvec,level,bb,dir,time);
00223 BBox wholebbox = base::GH().wholebbox();
00224 BBox inside = bb*coarsen(wholebbox,wholebbox.stepsize()/wholebbox.stepsize());
00225 if (!inside.empty())
00226 LBM().BCStandard(fvec,bb,LBMType::NoSlipWall,dir);
00227
00228 int scaling = LBMType::Physical;
00229 int Nx = fvec.extents(0), Ny = fvec.extents(1);
00230 int b = fvec.bottom();
00231 assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1));
00232 int mys = std::max(bb.lower(1)/bb.stepsize(1),fvec.lower(1)/fvec.stepsize(1));
00233 int mzs = std::max(bb.lower(2)/bb.stepsize(2),fvec.lower(2)/fvec.stepsize(2));
00234 int mxe = std::min(bb.upper(0)/bb.stepsize(0),fvec.upper(0)/fvec.stepsize(0));
00235 int mye = std::min(bb.upper(1)/bb.stepsize(1),fvec.upper(1)/fvec.stepsize(1));
00236 int mze = std::min(bb.upper(2)/bb.stepsize(2),fvec.upper(2)/fvec.stepsize(2));
00237 MicroType *f = (MicroType *)fvec.databuffer();
00238
00239 DataType rho=caux[0], a0=caux[1], a1=caux[2], a2=caux[3], elevation=caux[4], jr, hr;
00240
00241 MacroType q, q_;
00242 DataType profile, height, ldx;
00243
00244 if (scaling==LBMType::Physical) {
00245 rho /= LBM().DensityScale();
00246 a0 /= LBM().VelocityScale();
00247 a1 /= LBM().VelocityScale();
00248 a2 /= LBM().VelocityScale();
00249 elevation /= LBM().L0();
00250 }
00251
00252 if (bb.upper(0)+bb.stepsize(0)==wholebbox.lower(0)) {
00253
00254
00255 BBox wb = base::GH().wholebbox();
00256 DCoords wc;
00257 q_(0) = rho;
00258 q_(1) = a0;
00259 q_(2) = a1;
00260 q_(3) = a2;
00261
00262 int lc_indx[3];
00263
00264 lc_indx[0] = base::GH().worldCoords(fvec.lower(),fvec.stepsize())(0);
00265 lc_indx[1] = base::GH().worldCoords(fvec.lower(),fvec.stepsize())(1);
00266 ldx = DataType(fvec.stepsize(2))/DataType(wb.stepsize(2))*DataType(LBM().L0());
00267
00268 for (register int k=mzs; k<=mze; k++) {
00269 lc_indx[2] = k;
00270
00271 wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
00272 height = DataType(wc(2))*DataType(wb.stepsize(2));
00273 if ( height < elevation ) {
00274 profile = DataType(height*height)/DataType(elevation*elevation);
00275
00276 q(0) = q_(0);
00277 q(1) = q_(1)*profile;
00278 q(2) = q_(2)*profile;
00279 q(3) = q_(3)*profile;
00280 }
00281 else {
00282 q = q_;
00283 }
00284 for (register int j=mys; j<=mye; j++) {
00285 f[b+LBM().idx(mxe,j,k,Nx,Ny)] = LBM().Equilibrium(q);
00286 }
00287
00288
00289 }
00290 }
00291 }
00292 };
00293
00294 template <class LBMType, int dim>
00295 class LBMELCGFMBoundary : public LBMGFMBoundary<LBMType,dim> {
00296 typedef typename LBMType::MicroType MicroType;
00297 typedef typename MicroType::InternalDataType DataType;
00298 typedef LBMGFMBoundary<LBMType,dim> base;
00299 public:
00300 typedef typename base::vec_grid_data_type vec_grid_data_type;
00301 typedef typename base::grid_data_type grid_data_type;
00302 typedef typename base::point_type point_type;
00303 typedef AMRELCGFMSolver<MicroType,FixupType,FlagType,dim> amr_elc_solver;
00304
00305 LBMELCGFMBoundary(LBMType &scheme, amr_elc_solver& solver) : base(scheme), _solver(solver) {
00306 base::SetNAux(base::Dim());
00307 }
00308
00309 virtual void SetBndryAux(vec_grid_data_type& gdu, grid_data_type& gdphi,
00310 const MicroType* u, DataType* aux, const int& Level,
00311 double t, const int& nc, const int* idx,
00312 const point_type* xc, const DataType* distance,
00313 const point_type* normal) {
00314 _solver.ExtractBoundaryVelocities(gdu.bbox(),nc,xc,reinterpret_cast<point_type *>(aux));
00315 }
00316
00317 protected:
00318 amr_elc_solver& _solver;
00319 };
00320
00321 class SolverObjects {
00322 public:
00323 SolverObjects() {
00324 _LBMSpecific = new LBMSpecific();
00325 _IntegratorSpecific = new IntegratorSpecific(*_LBMSpecific);
00326 _InitialConditionSpecific = new InitialConditionSpecific(*_LBMSpecific);
00327 _BoundaryConditionsSpecific = new BoundaryConditionsSpecific(*_LBMSpecific);
00328 }
00329
00330 ~SolverObjects() {
00331 delete _BoundaryConditionsSpecific;
00332 delete _InitialConditionSpecific;
00333 delete _IntegratorSpecific;
00334 delete _LBMSpecific;
00335 }
00336
00337 protected:
00338 LBMSpecific *_LBMSpecific;
00339 IntegratorSpecific *_IntegratorSpecific;
00340 InitialConditionSpecific *_InitialConditionSpecific;
00341 BoundaryConditionsSpecific *_BoundaryConditionsSpecific;
00342 };
00343
00344 class FluidSolverSpecific : protected SolverObjects,
00345 public AMRELCGFMSolver<MicroType,FixupType,FlagType,DIM> {
00346 typedef MicroType::InternalDataType DataType;
00347 typedef AMRELCGFMSolver<MicroType,FixupType,FlagType,DIM> base;
00348 typedef AMRInterpolation<MicroType,DIM> interpolation_type;
00349 typedef interpolation_type::point_type point_type;
00350 typedef SchemeGFMFileOutput<LBMType,FixupType,FlagType,DIM> output_type;
00351 typedef ads::FixedArray<10,DataType> PVCType;
00352
00353 typedef AMRStatistics<MicroType,interpolation_type,output_type,DIM> stat_type;
00354 public:
00355 FluidSolverSpecific() : SolverObjects(), base(*SolverObjects::_IntegratorSpecific,
00356 *SolverObjects::_InitialConditionSpecific,
00357 *SolverObjects::_BoundaryConditionsSpecific) {
00358 SetLevelTransfer(new LBMF77LevelTransfer<LBMType,DIM>(SolverObjects::_IntegratorSpecific->Scheme(), f_prolong, f_restrict));
00359 SetFileOutput(new output_type(*this,SolverObjects::_IntegratorSpecific->Scheme()));
00360 SetFixup(new FixupSpecific(SolverObjects::_IntegratorSpecific->Scheme()));
00361 SetFlagging(new FlaggingSpecific(*this,SolverObjects::_IntegratorSpecific->Scheme()));
00362 AddGFM(new GhostFluidMethod<MicroType,DIM>(new LBMELCGFMBoundary<LBMType,DIM>(SolverObjects::_IntegratorSpecific->LBM(),*this),
00363 new CPTLevelSet<DataType,DIM>()));
00364
00365 SetCoupleGFM(0);
00366
00367 _Stats = new stat_type(_Interpolation, (output_type*)_FileOutput);
00368 }
00369
00370 ~FluidSolverSpecific() {
00371 DeleteGFM(_GFM[0]);
00372 delete _LevelTransfer;
00373 delete _Flagging;
00374 delete _Fixup;
00375 delete _FileOutput;
00376 delete _BoundaryConditionsSpecific;
00377 _BoundaryConditionsSpecific = (BoundaryConditionsSpecific *) 0;
00378 delete _InitialConditionSpecific;
00379 _InitialConditionSpecific = (InitialConditionSpecific *) 0;
00380 delete _IntegratorSpecific;
00381 _IntegratorSpecific = (IntegratorSpecific *) 0;
00382 delete _Interpolation;
00383 _Interpolation = (AMRELCGFMSolver<MicroType,FixupType,FlagType,DIM>::interpolation_type *) 0;
00384 delete _Stats;
00385 _Stats = (stat_type *) 0;
00386 }
00387
00388 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00389 base::register_at(Ctrl,prefix);
00390 if (_Stats) _Stats->register_at(base::LocCtrl, "");
00391 RegisterAt(base::LocCtrl,"track_every",track_every);
00392 }
00393 virtual void register_at(ControlDevice& Ctrl) {
00394 base::register_at(Ctrl);
00395 }
00396
00397 virtual void SendBoundaryData() {
00398 START_WATCH
00399 for (register int l=0; l<=FineLevel(base::GH()); l++) {
00400 int Time = CurrentTime(base::GH(),l);
00401 ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time,
00402 l, base::Dim()+4, base::t[l]);
00403 if (CurrentTime(base::GH(),base::CouplingLevel) != Time)
00404 ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time+TimeStep(base::U(),l),
00405 l, base::Dim()+4, base::t[l]+base::dt[l]);
00406 }
00407 END_WATCH(FLUID_CPL_PRESSURE_CALCULATE)
00408 base::SendBoundaryData();
00409 }
00410
00411 virtual void SetupData() {
00412 base::SetupData();
00413 base::NAMRTimeSteps = 1;
00414 base::AdaptBndTimeInterpolate = 0;
00415 base::Step[0].LastTime = 1.e37;
00416 base::Step[0].VariableTimeStepping = -1;
00417 base::Step[0].dtv[0] = ((LBMIntegrator<LBMType,DIM> &)Integrator_()).LBM().TimeScale();
00418
00419 _Interpolation->SetupData(base::PGH(), base::NGhosts());
00420 _Stats->Setup(base::PGH(), base::NGhosts(), base::shape, base::geom);
00421 }
00422
00423 virtual void Advance(double& t, double& dt) {
00424 base::Advance(t,dt);
00425 _Stats->Evaluate(base::t, base::U(), base::Work());
00426
00427 }
00428
00429 virtual void Initialize_(const double& dt_start) {
00430 base::Initialize_(dt_start);
00431 _Stats->Evaluate(base::t, base::U(), base::Work());
00432 }
00433
00434 double modulus(double a, double b) const {
00435 int result = (int) std::floor( a / b );
00436 return a - static_cast<double>( result ) * b;
00437 }
00438
00439 bool inBox (DCoords min, DCoords max, DCoords wc) const {
00440 if ( min(0) <= wc(0) && wc(0) <= max(0)) {
00441
00442 if ( min(1) <= wc(1) && wc(1) <= max(1)) {
00443
00444 if ( min(2) <= wc(2) && wc(2) <= max(2)) {
00445 return 1;
00446 }
00447 }
00448 }
00449 return 0;
00450 }
00451
00452 protected:
00453 int track_every;
00454 stat_type* _Stats;
00455 };