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 #define OWN_INITIALCONDITION
00016 #define OWN_BOUNDARYCONDITION
00017 #include "LBMStdProblem.h"
00018 #include "F77Interfaces/F77GFMBoundary.h"
00019 #include "AMRELCGFMSolver.h"
00020 #include "LBMGFMBoundary.h"
00021 #include "Interfaces/SchemeGFMFileOutput.h"
00022 #include "AMRInterpolation.h"
00023 #define NVARS 5
00024 #include "AMRStatistics.h"
00025
00026 class LBMSpecific : public LBMType {
00027 typedef LBMType base;
00028 public:
00029 LBMSpecific() : base(), omega(1.), u_p_avg(1.), l0_p(-1.) {}
00030
00031 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00032 base::register_at(Ctrl,prefix);
00033 RegisterAt(base::LocCtrl,"Omega",omega);
00034 RegisterAt(base::LocCtrl,"u_p_avg",u_p_avg);
00035 RegisterAt(base::LocCtrl,"l0_p",l0_p);
00036 }
00037
00038 virtual void SetupData(GridHierarchy* gh, const int& ghosts) {
00039 base::SetupData(gh,ghosts);
00040
00041 const double* bndry = base::GH().wholebndry();
00042 const int nbndry = base::GH().nbndry();
00043 if (l0_p <= 0. ) {
00044 l0_p=bndry[nbndry*(1+2*1)]-bndry[nbndry*(0+2*1)];
00045 if (nbndry>1) l0_p=bndry[nbndry*(1+2*1)]-bndry[1+nbndry*(1+2*1)];
00046 }
00047 DataType omega = base::Omega(base::TimeScale());
00048 double Re_p=u_p_avg*l0_p/base::GasViscosity();
00049 double Re_lbm=(u_p_avg/base::VelocityScale())*(l0_p/base::LengthScale())/base::LatticeViscosity(omega);
00050 ( comm_service::log() << "l0_p=" << l0_p << " u_p_avg=" << u_p_avg
00051 << " Re_p=" << Re_p << " Re_lbm=" << Re_lbm
00052 << "\n omega=" << omega
00053 << " SpeedUp=" << base::SpeedUp()
00054 << " u_avg_LBM=" << u_p_avg/base::VelocityScale()
00055 << " SpeedRatio=" << u_p_avg/base::VelocityScale()*std::sqrt(DataType(3.))
00056 << " base::TO()=" << base::T0()
00057 << "\n Re_p-Re_lbm=" << Re_p-Re_lbm
00058 << std::endl ).flush();
00059 if ( std::fabs(Re_p-Re_lbm)>DBL_EPSILON*LB_FACTOR*LB_FACTOR*LB_FACTOR ) {
00060 std::cerr << "Physical scaling requires fabs(Re_p-Re_lbm)<DBL_EPSILON*LB_FACTOR*LB_FACTOR*LB_FACTOR !" << std::endl;
00061 std::cerr << "Re_p " << Re_p << " Re_lbm " << Re_lbm << " std::fabs(Re_p-Re_lbm) " << std::fabs(Re_p-Re_lbm) << std::endl;
00062 std::exit(-1);
00063 }
00064 if (omega<0.5 || omega>2.) {
00065 std::cerr << "0.5 <= omega <= 2 required." << std::endl;
00066 std::exit(-1);
00067 }
00068 base::SetGas(base::GasDensity(),base::GasViscosity(),base::GasSpeedofSound());
00069 ( comm_service::log() << "FluidProblem SetupData() D3Q19: Gas_rho=" << base::GasDensity() << " Gas_Cs=" << base::GasSpeedofSound()
00070 << " Gas_nu=" << base::GasViscosity() << std::endl ).flush();
00071
00072 }
00073
00074 private:
00075 DataType omega, u_p_avg, l0_p;
00076 };
00077
00078 class InitialConditionSpecific : public LBMInitialCondition<LBMType,DIM> {
00079 typedef LBMInitialCondition<LBMType,DIM> base;
00080 public:
00081 InitialConditionSpecific(LBMType &lbm) : base(lbm) {}
00082 typedef base::MicroType MicroType;
00083 typedef base::MacroType MacroType;
00084
00085 virtual void SetGrid(vec_grid_data_type& fvec, grid_data_type& workvec, const int& level) {
00086 MacroType q, q_;
00087 DataType profile, height;
00088 BBox bb = fvec.bbox();
00089 DCoords dx = base::GH().worldStep(fvec.stepsize());
00090
00091 q_(0) = base::aux[0];
00092 q_(1) = base::aux[1];
00093 q_(2) = base::aux[2];
00094 q_(3) = base::aux[3];
00095 DataType elevation = base::aux[4];
00096
00097 if (base::Scaling==LBMType::Physical) {
00098 q_(0) /= LBM().DensityScale();
00099 q_(1) /= LBM().VelocityScale();
00100 q_(2) /= LBM().VelocityScale();
00101 q_(3) /= LBM().VelocityScale();
00102 }
00103
00104 MicroType *f = (MicroType *)fvec.databuffer();
00105
00106 int Nx = fvec.extents(0), Ny = fvec.extents(1);
00107 int b = fvec.bottom();
00108 assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1));
00109 int mxs = std::max(bb.lower(0)/bb.stepsize(0),fvec.lower(0)/fvec.stepsize(0));
00110 int mys = std::max(bb.lower(1)/bb.stepsize(1),fvec.lower(1)/fvec.stepsize(1));
00111 int mzs = std::max(bb.lower(2)/bb.stepsize(2),fvec.lower(2)/fvec.stepsize(2));
00112 int mxe = std::min(bb.upper(0)/bb.stepsize(0),fvec.upper(0)/fvec.stepsize(0));
00113 int mye = std::min(bb.upper(1)/bb.stepsize(1),fvec.upper(1)/fvec.stepsize(1));
00114 int mze = std::min(bb.upper(2)/bb.stepsize(2),fvec.upper(2)/fvec.stepsize(2));
00115 int lc_indx[3] = { mxs*fvec.stepsize(0), mys*fvec.stepsize(1), 0 };
00116
00117 for (register int k=mzs; k<=mze; k++) {
00118 lc_indx[2] = k*fvec.stepsize(2);
00119 DCoords wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
00120 height = wc(2)+dx(2)/2.;
00121
00122 if ( height < elevation ) {
00123 profile = (height*height)/(elevation*elevation);
00124 q(0) = q_(0);
00125 q(1) = q_(1)*profile;
00126 q(2) = q_(2)*profile;
00127 q(3) = q_(3)*profile;
00128 }
00129 else
00130 q = q_;
00131
00132 for (register int j=mys; j<=mye; j++)
00133 for (register int i=mxs; i<=mxe; i++)
00134 f[b+LBM().idx(i,j,k,Nx,Ny)] = LBM().Equilibrium(q);
00135 }
00136 }
00137 };
00138
00139 class BoundaryConditionsSpecific : public LBMBoundaryConditions<LBMType,DIM> {
00140 typedef LBMBoundaryConditions<LBMType,DIM> base;
00141 public:
00142 typedef base::MicroType MicroType;
00143 typedef base::MacroType MacroType;
00144 BoundaryConditionsSpecific(LBMType &lbm) : base(lbm) {}
00145
00146 virtual void SetupData(GridHierarchy* gh, const int& ghosts) {
00147 base::SetupData(gh,ghosts);
00148 const double* bndry = base::GH().wholebndry();
00149 const int nbndry = base::GH().nbndry();
00150 DataType l0_p=bndry[nbndry*(1+2*1)]-bndry[nbndry*(0+2*1)];
00151 if (nbndry>1) l0_p=bndry[nbndry*(1+2*1)]-bndry[1+nbndry*(1+2*1)];
00152 DataType omega = LBM().Omega(LBM().TimeScale());
00153 DataType u_p_avg=Side(0).aux[0];
00154 double Re_p=u_p_avg*l0_p/LBM().GasViscosity();
00155 double Re_lbm=(u_p_avg/LBM().VelocityScale())*(l0_p/LBM().LengthScale())/LBM().LatticeViscosity(omega);
00156 ( comm_service::log() << "l0_p=" << l0_p << " u_p_avg=" << u_p_avg
00157 << " Re_p=" << Re_p << " Re_lbm=" << Re_lbm
00158 << "\n omega=" << omega
00159 << " SpeedUp=" << LBM().SpeedUp()
00160 << " u_avg_LBM=" << u_p_avg/LBM().VelocityScale()
00161 << " SpeedRatio=" << u_p_avg/LBM().VelocityScale()*std::sqrt(DataType(3.))
00162 << " LBM().TO()=" << LBM().T0()
00163 << "\n Re_p-Re_lbm=" << Re_p-Re_lbm
00164 << std::endl ).flush();
00165 assert(std::fabs(Re_p-Re_lbm)<DBL_EPSILON*LB_FACTOR*LB_FACTOR*LB_FACTOR);
00166 }
00167
00168 virtual void SetBndry(vec_grid_data_type &fvec, const int& level, const BBox &bb,
00169 const int &dir, const double& time) {
00170
00171 base::SetBndry(fvec,level,bb,dir,time);
00172 BBox wholebbox = base::GH().wholebbox();
00173 BBox inside = bb*coarsen(wholebbox,wholebbox.stepsize()/wholebbox.stepsize());
00174 if (!inside.empty())
00175 LBM().BCStandard(fvec,bb,LBMType::NoSlipWall,dir);
00176
00177 if (dir==LBMType::Left && bb.upper(0)+bb.stepsize(0)==wholebbox.lower(0) &&
00178 Side(0).Type==7) {
00179
00180 int Nx = fvec.extents(0), Ny = fvec.extents(1);
00181 int b = fvec.bottom();
00182 assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1)
00183 && fvec.stepsize(2)==bb.stepsize(2));
00184 int mxs = std::max(bb.lower(0)/bb.stepsize(0),fvec.lower(0)/fvec.stepsize(0));
00185 int mys = std::max(bb.lower(1)/bb.stepsize(1),fvec.lower(1)/fvec.stepsize(1));
00186 int mzs = std::max(bb.lower(2)/bb.stepsize(2),fvec.lower(2)/fvec.stepsize(2));
00187 int mxe = std::min(bb.upper(0)/bb.stepsize(0),fvec.upper(0)/fvec.stepsize(0));
00188 int mye = std::min(bb.upper(1)/bb.stepsize(1),fvec.upper(1)/fvec.stepsize(1));
00189 int mze = std::min(bb.upper(2)/bb.stepsize(2),fvec.upper(2)/fvec.stepsize(2));
00190 MicroType *f = (MicroType *)fvec.databuffer();
00191
00192 DataType rho=Side(0).aux[0];
00193 DataType a0=Side(0).aux[1];
00194 DataType a1=Side(0).aux[2];
00195 DataType a2=Side(0).aux[3];
00196 DataType elevation=Side(0).aux[4];
00197
00198 MacroType q, q_;
00199 DataType profile, height;
00200 DCoords dx = base::GH().worldStep(fvec.stepsize());
00201 int lc_indx[3] = { mxs*fvec.stepsize(0), mys*fvec.stepsize(1), 0 };
00202
00203 if (Side(0).Scaling==LBMType::Physical) {
00204 rho /= LBM().DensityScale();
00205 a0 /= LBM().VelocityScale();
00206 a1 /= LBM().VelocityScale();
00207 a2 /= LBM().VelocityScale();
00208 }
00209
00210 q_(0) = rho;
00211 q_(1) = a0;
00212 q_(2) = a1;
00213 q_(3) = a2;
00214
00215 for (register int k=mzs; k<=mze; k++) {
00216 lc_indx[2] = k*fvec.stepsize(2);
00217 DCoords wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
00218 height = wc(2)+dx(2)/2.;
00219
00220 if ( height < elevation ) {
00221 profile = DataType(height*height)/DataType(elevation*elevation);
00222 q(0) = q_(0);
00223 q(1) = q_(1)*profile;
00224 q(2) = q_(2)*profile;
00225 q(3) = q_(3)*profile;
00226 }
00227 else
00228 q = q_;
00229
00230 for (register int j=mys; j<=mye; j++)
00231 f[b+LBM().idx(mxe,j,k,Nx,Ny)] = LBM().Equilibrium(q);
00232 }
00233 }
00234 }
00235 };
00236
00237
00238 template <class LBMType, int dim>
00239 class LBMELCGFMBoundary : public LBMGFMBoundary<LBMType,dim> {
00240 typedef typename LBMType::MicroType MicroType;
00241 typedef typename MicroType::InternalDataType DataType;
00242 typedef LBMGFMBoundary<LBMType,dim> base;
00243 public:
00244 typedef typename base::vec_grid_data_type vec_grid_data_type;
00245 typedef typename base::grid_data_type grid_data_type;
00246 typedef typename base::point_type point_type;
00247 typedef AMRELCGFMSolver<MicroType,FixupType,FlagType,dim> amr_elc_solver;
00248
00249 LBMELCGFMBoundary(LBMType &scheme, amr_elc_solver& solver) : base(scheme), _solver(solver) {
00250 base::SetNAux(base::Dim());
00251 }
00252
00253 virtual void SetBndryAux(vec_grid_data_type& gdu, grid_data_type& gdphi,
00254 const MicroType* u, DataType* aux, const int& Level,
00255 double t, const int& nc, const int* idx,
00256 const point_type* xc, const DataType* distance,
00257 const point_type* normal) {
00258 _solver.ExtractBoundaryVelocities(gdu.bbox(),nc,xc,reinterpret_cast<point_type *>(aux));
00259 }
00260
00261 protected:
00262 amr_elc_solver& _solver;
00263 };
00264
00265 class SolverObjects {
00266 public:
00267 SolverObjects() {
00268 _LBMSpecific = new LBMSpecific();
00269 _IntegratorSpecific = new IntegratorSpecific(*_LBMSpecific);
00270 _InitialConditionSpecific = new InitialConditionSpecific(*_LBMSpecific);
00271 _BoundaryConditionsSpecific = new BoundaryConditionsSpecific(*_LBMSpecific);
00272 }
00273
00274 ~SolverObjects() {
00275 delete _BoundaryConditionsSpecific;
00276 delete _InitialConditionSpecific;
00277 delete _IntegratorSpecific;
00278 delete _LBMSpecific;
00279 }
00280
00281 protected:
00282 LBMSpecific *_LBMSpecific;
00283 IntegratorSpecific *_IntegratorSpecific;
00284 InitialConditionSpecific *_InitialConditionSpecific;
00285 BoundaryConditionsSpecific *_BoundaryConditionsSpecific;
00286 };
00287
00288 class FluidSolverSpecific : protected SolverObjects,
00289 public AMRELCGFMSolver<MicroType,FixupType,FlagType,DIM> {
00290 typedef MicroType::InternalDataType DataType;
00291 typedef AMRELCGFMSolver<MicroType,FixupType,FlagType,DIM> base;
00292 typedef SchemeGFMFileOutput<LBMType,FixupType,FlagType,DIM> output_type;
00293 typedef AMRInterpolation<MicroType,DIM> interpolation_type;
00294 typedef interpolation_type::point_type point_type;
00295 typedef AMRStatistics<MicroType,interpolation_type,output_type,DIM> stat_type;
00296 public:
00297 FluidSolverSpecific() : SolverObjects(), base(*SolverObjects::_IntegratorSpecific,
00298 *SolverObjects::_InitialConditionSpecific,
00299 *SolverObjects::_BoundaryConditionsSpecific) {
00300 SetLevelTransfer(new LBMF77LevelTransfer<LBMType,DIM>(SolverObjects::_IntegratorSpecific->Scheme(), f_prolong, f_restrict));
00301 SetFileOutput(new output_type(*this,SolverObjects::_IntegratorSpecific->Scheme()));
00302
00303 SetFixup(new FixupSpecific(SolverObjects::_IntegratorSpecific->Scheme()));
00304 SetFlagging(new FlaggingSpecific(*this,SolverObjects::_IntegratorSpecific->Scheme()));
00305 AddGFM(new GhostFluidMethod<MicroType,DIM>(new LBMELCGFMBoundary<LBMType,DIM>(SolverObjects::_IntegratorSpecific->LBM(),*this),
00306 new CPTLevelSet<DataType,DIM>()));
00307 SetCoupleGFM(0);
00308 _Interpolation = new interpolation_type();
00309 _Stats = new stat_type(_Interpolation, (output_type*)_FileOutput);
00310 }
00311
00312 ~FluidSolverSpecific() {
00313 DeleteGFM(_GFM[0]);
00314 delete _LevelTransfer;
00315 delete _Flagging;
00316 delete _Fixup;
00317 delete _FileOutput;
00318 delete _Interpolation;
00319 delete _Stats;
00320 }
00321
00322 virtual void register_at(ControlDevice& Ctrl) { base::register_at(Ctrl); }
00323 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00324 base::register_at(Ctrl, prefix);
00325 if (_Stats) _Stats->register_at(base::LocCtrl, "");
00326 }
00327
00328 virtual void SetupData() {
00329 ( comm_service::log() << " fluid setup begun" << std::endl ).flush();
00330 base::SetupData();
00331 base::AdaptBndTimeInterpolate = 0;
00332 base::NAMRTimeSteps = 1;
00333 base::Step[0].LastTime = 1.e37;
00334 base::Step[0].VariableTimeStepping = -1;
00335 base::Step[0].dtv[0] = ((LBMIntegrator<LBMType,DIM> &)Integrator_()).LBM().TimeScale();
00336 _Interpolation->SetupData(base::PGH(), base::NGhosts());
00337 _Stats->Setup(base::PGH(), base::NGhosts(), base::shape, base::geom);
00338 }
00339
00340 virtual void SendBoundaryData() {
00341 START_WATCH
00342 for (register int l=0; l<=FineLevel(base::GH()); l++) {
00343 int Time = CurrentTime(base::GH(),l);
00344 ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time,
00345 l, base::Dim()+4, base::t[l]);
00346 if (CurrentTime(base::GH(),base::CouplingLevel) != Time)
00347 ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time+TimeStep(base::U(),l),
00348 l, base::Dim()+4, base::t[l]+base::dt[l]);
00349 }
00350 END_WATCH(FLUID_CPL_PRESSURE_CALCULATE)
00351 base::SendBoundaryData();
00352 }
00353
00354 virtual void WaitReceiveBoundaryData() {
00355 base::WaitReceiveBoundaryData();
00356 GFM(0).LevelSet().SetStationary(0);
00357 }
00358
00359 virtual void RecomposeGridHierarchy(const int Time, const int Level,
00360 bool ShadowAllowed, bool DoFixup,
00361 bool RecomposeBaseLev, bool RecomposeHighLev) {
00362 base::RecomposeGridHierarchy(Time,Level,ShadowAllowed,DoFixup,RecomposeBaseLev,RecomposeHighLev);
00363 GFM(0).LevelSet().SetStationary(1);
00364 }
00365
00366 virtual void Advance(double& t, double& dt) {
00367 base::Advance(t, dt);
00368 _Stats->Evaluate(base::t, base::U(), base::Work());
00369 }
00370
00371 virtual void Initialize_(const double& dt_start) {
00372 base::Initialize_(dt_start);
00373 _Stats->Evaluate(base::t, base::U(), base::Work());
00374 }
00375
00376 protected:
00377 interpolation_type* _Interpolation;
00378 stat_type* _Stats;
00379 };