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