00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #ifndef AMROC_PROBLEM_H
00013 #define AMROC_PROBLEM_H
00014
00015 #define DIM 2
00016
00017 #include "LBMProblem.h"
00018 #include "LBMD2Q9Thermal.h"
00019
00020 typedef double DataType;
00021 typedef LBMD2Q9Thermal<DataType> LBMType;
00022
00023 #define OWN_LBMSCHEME
00024 #define OWN_AMRSOLVER
00025 #include "LBMStdProblem.h"
00026 #include "AMRGFMSolver.h"
00027 #include "AMRGFMInterpolation.h"
00028 #include "Interfaces/SchemeGFMFileOutput.h"
00029
00030 DataType rd, xs[2], TempS, vr;
00031 DataType pi=4.0*std::atan(1.0);
00032
00033 class LBMSpecific : public LBMType {
00034 typedef LBMType base;
00035 public:
00036 LBMSpecific() : base(), Tmp(1.), deltaTp(1.), gp(1.), betap(1.), Diffp(1.),GasRho(1.), Gasnu(1.), Gasp0(1.), Ra(1.) {}
00037
00038 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00039 base::register_at(Ctrl,prefix);
00040 RegisterAt(base::LocCtrl,"Tm",Tmp);
00041 RegisterAt(base::LocCtrl,"deltaT",deltaTp);
00042 RegisterAt(base::LocCtrl,"g",gp);
00043 RegisterAt(base::LocCtrl,"beta",betap);
00044 RegisterAt(base::LocCtrl,"Diff",Diffp);
00045 RegisterAt(base::LocCtrl,"Gas_rho",GasRho);
00046 RegisterAt(base::LocCtrl,"Gas_nu",Gasnu);
00047 RegisterAt(base::LocCtrl,"Gas_p0",Gasp0);
00048 RegisterAt(base::LocCtrl,"Ra",Ra);
00049 }
00050
00051 virtual void SetupData(GridHierarchy* gh, const int& ghosts) {
00052 base::SetupData(gh,ghosts);
00053
00054 BBox wb = base::GH().wholebbox();
00055 DCoords lc = base::GH().worldCoords(wb.lower(),wb.stepsize());
00056 DCoords uc = base::GH().worldCoords(wb.upper()+wb.stepsize(),wb.stepsize());
00057
00058
00059
00060 DataType Hp=uc(1)-lc(1);
00061
00062
00063
00064 DataType viscp = Gasnu;
00065 DataType Pr = viscp/Diffp;
00066
00067 DataType tScale = base::TimeScale();
00068 DataType lScale = base::LengthScale();
00069 DataType THp = Tmp+deltaTp/DataType(2.);
00070 DataType TCp = Tmp-deltaTp/DataType(2.);
00071 DataType Ma = std::sqrt(DataType(3.)*gp*betap*deltaTp*Hp);
00072 DataType MaS = std::sqrt(DataType(3.)*Ra/Pr)*viscp*base::SpeedUp()/Hp;
00073
00074
00075 DataType GasCsp = base::LatticeSpeedOfSound();
00076
00077 base::SetGas(Gasp0,GasRho,viscp,GasCsp);
00078 base::SetThermalGas(TCp,THp,Diffp,gp,betap);
00079
00080 std::cout << "D2Q19Thermal: Gas_rho=" << base::GasDensity() << " Gas_Cs=" << base::GasSpeedofSound()
00081 << " Gas_nu=" << base::GasViscosity() << " Gas_nut=" << base::GasViscosityT() << std::endl;
00082 std::cout << "\nControl Parameters : " << std::endl;
00083 std::cout << "Ra = " << gp*betap*Hp*Hp*Hp/viscp/Diffp << " Pr = " << Pr << std::endl;
00084 std::cout << "\nTemperature : " << std::endl;
00085 std::cout << "TH = " << THp << " TC = " << TCp << std::endl;
00086 std::cout << "\nFluid values : " << std::endl;
00087 std::cout << "visc = " << viscp << " Diff = " << Diffp << " g = " << gp << " betap = " << betap << std::endl;
00088 std::cout << "Ma = " << Ma << " MaS = " << MaS << std::endl;
00089 std::cout << "\nGeometry : " << std::endl;
00090 std::cout << "Channel Width = " << Hp << std::endl;
00091 std::cout << "\nLattice Units : " << std::endl;
00092 std::cout << "(Amroc) Omega = " << base::Omega(base::TimeScale()) << " (Amroc) OmegaT = " << base::OmegaT(base::TimeScale()) << std::endl;
00093 std::cout << "(Amroc) Gas_nu = " << base::LatticeViscosity(base::Omega(tScale)) << " (Amroc) Gas_nuT = " << base::LatticeViscosityT(base::OmegaT(tScale)) << std::endl;
00094
00095
00096
00097 std::cout << "\nScaling Factors: " << std::endl;
00098 std::cout << "LengthScale = " << base::LengthScale() << std::endl;
00099 std::cout << "TimeScale = " << base::TimeScale() << std::endl;
00100 std::cout << "VelocityScale = " << base::VelocityScale() << std::endl;
00101 std::cout << "SpeedUp = " << base::SpeedUp() << std::endl;
00102
00103 }
00104
00105 protected:
00106 DataType Tmp, deltaTp, gp, betap, Diffp;
00107 DataType GasRho, Gasnu, Gasp0, Ra;
00108
00109 };
00110
00111 class LBMExactSolutionSpecific : public LBMExactSolution<LBMType,DIM> {
00112 typedef LBMExactSolution<LBMType,DIM> base;
00113 typedef AMRGFMSolver<MicroType,FixupType,FlagType,DIM> solver_type;
00114 public:
00115 typedef base::vec_grid_fct_type vec_grid_fct_type;
00116 typedef base::vec_grid_data_type vec_grid_data_type;
00117 typedef LBMType::MicroType MicroType;
00118 typedef LBMType::MacroType MacroType;
00119
00120 LBMExactSolutionSpecific(LBMType &lbm) : base(lbm) {
00121
00122 u0=0.; TH0=0.; deltaT0=0.; L0=0.; Re0=0.; Pr0=0.;
00123 }
00124
00125 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00126 base::register_at(Ctrl,prefix);
00127 RegisterAt(base::LocCtrl,"u0",u0);
00128 RegisterAt(base::LocCtrl,"TH0",TH0);
00129 RegisterAt(base::LocCtrl,"deltaT0",deltaT0);
00130 RegisterAt(base::LocCtrl,"L0",L0);
00131 RegisterAt(base::LocCtrl,"Re0",Re0);
00132 RegisterAt(base::LocCtrl,"Pr0",Pr0);
00133 }
00134
00135 virtual void SetGrid(vec_grid_data_type& fvec, grid_data_type& workvec, const int& level, const double& t) {
00136
00137 int Nx = fvec.extents(0), Ny = fvec.extents(1);
00138 int mxs = base::NGhosts(), mys = base::NGhosts();
00139 int mxe = Nx-base::NGhosts()-1, mye = Ny-base::NGhosts()-1;
00140
00141 MicroType *f = (base::MicroType *)fvec.databuffer();
00142
00143
00144
00145 DCoords lbcorner = base::GH().worldCoords(fvec.lower(), fvec.stepsize());
00146 DCoords dx = base::GH().worldStep(fvec.stepsize());
00147
00148 for (register int j=mys; j<=mye; j++){
00149 for (register int i=mxs; i<=mxe; i++) {
00150
00151
00152 double xc = (i+0.5)*dx(0)+lbcorner(0);
00153 double yc = (j+0.5)*dx(1)+lbcorner(1);
00154
00155
00156 MacroType q = LBM().MacroVariables(f[LBM().idx(i,j,Nx)]);
00157 MacroType qr = q;
00158
00159
00160 qr(1) = u0*LBM().SpeedUp()*((std::exp(Re0*yc/L0)-1.)/(std::exp(Re0)-1.));
00161 qr(3) = TH0-deltaT0*((std::exp(Re0*Pr0*yc/L0)-1.)/(std::exp(Re0*Pr0)-1.));
00162
00163
00164 for (register int l=0; l<LBM().NMacroVar(); l++)
00165 f[LBM().idx(i,j,Nx)](l) += std::abs(q(l)-qr(l));
00166 }
00167 }
00168
00169 }
00170
00171 inline LBMType& LBM() { return base::Scheme(); }
00172 inline const LBMType& LBM() const { return base::Scheme(); }
00173
00174 protected:
00175 double u0, TH0, deltaT0, L0, Re0, Pr0;
00176 };
00177
00178
00179 class SolverSpecific :
00180 public AMRGFMSolver<MicroType,FixupType,FlagType,DIM> {
00181 typedef AMRGFMSolver<MicroType,FixupType,FlagType,DIM> base;
00182 typedef AMRGFMInterpolation<MicroType,FixupType,FlagType,DIM> interpolation_type;
00183 typedef MicroType::InternalDataType DataType;
00184 typedef SchemeGFMFileOutput<LBMType,FixupType,FlagType,DIM> output_type;
00185 typedef interpolation_type::point_type point_type;
00186 public:
00187 SolverSpecific(IntegratorSpecific& integ,
00188 base::initial_condition_type& init,
00189 base::boundary_conditions_type& bc) : base(integ, init, bc) {
00190 SetLevelTransfer(new LBMF77LevelTransfer<LBMType,DIM>(integ.Scheme(), f_prolong, f_restrict));
00191 SetFileOutput(new SchemeGFMFileOutput<LBMType,FixupType,FlagType,DIM>(*this,integ.Scheme()));
00192 SetFixup(new FixupSpecific(integ.Scheme()));
00193 SetFlagging(new FlaggingSpecific(*this,integ.Scheme()));
00194 _Interpolation = new interpolation_type(*this);
00195 SetExactSolution(new LBMExactSolutionSpecific(integ.Scheme()));
00196
00197 u0=0.; TH0=0.; deltaT0=0.; L0=0.; Re0=0.; Pr0=0.; IntegrateEvery=10.;
00198
00199 }
00200
00201 ~SolverSpecific() {
00202 DeleteGFM(_GFM[0]);
00203 delete _LevelTransfer;
00204 delete _Flagging;
00205 delete _Fixup;
00206 delete _FileOutput;
00207 delete _Interpolation;
00208 delete _ExactSolution;
00209
00210 }
00211
00212 virtual void SetupData() {
00213 base::SetupData();
00214 base::NAMRTimeSteps = 1;
00215 base::AdaptBndTimeInterpolate = 0;
00216 base::Step[0].LastTime = 1.e37;
00217 base::Step[0].VariableTimeStepping = -1;
00218 base::Step[0].dtv[0] = ((LBMIntegrator<LBMType,DIM> &)Integrator_()).LBM().TimeScale();
00219 _Interpolation->SetupData(base::PGH(), base::NGhosts());
00220
00221 const double* bndry = base::GH().wholebndry();
00222 const int nbndry = base::GH().nbndry();
00223
00224 LBMType &lbm = (LBMType &)((LBMIntegrator<LBMType,DIM> &)Integrator_()).LBM();
00225 u0 = ((LBMBoundaryConditions<LBMType,DIM> &)BoundaryConditions_()).Side(3).aux[5];
00226 double v0 = ((LBMBoundaryConditions<LBMType,DIM> &)BoundaryConditions_()).Side(2).aux[6];
00227 double TC0 = ((LBMBoundaryConditions<LBMType,DIM> &)BoundaryConditions_()).Side(3).aux[7];
00228 TH0 = ((LBMBoundaryConditions<LBMType,DIM> &)BoundaryConditions_()).Side(2).aux[7];
00229 deltaT0 = TH0-TC0;
00230 L0 = bndry[nbndry*(1+2*1)]-bndry[nbndry*(0+2*1)];
00231 double omega0 = lbm.Omega(lbm.TimeScale());
00232 Re0 = v0*L0/lbm.GasViscosity();
00233 double Re_0LBM = (v0/lbm.VelocityScale())*(L0/lbm.LengthScale())/lbm.LatticeViscosity(omega0);
00234 Pr0 = lbm.GasViscosity()/lbm.GasViscosityT();
00235
00236 std::cout << "Debug: u0 = " << u0 << " TC0 = " << TC0 << " TH0 = " << TH0 << " deltaT0 = " << deltaT0 << " L0 = " << L0 << " Re0 = " << Re0 << " ReLBM0 = " << Re_0LBM << " Pr0 = " << Pr0 <<std::endl;
00237 }
00238
00239 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00240 base::register_at(Ctrl,prefix);
00241 RegisterAt(base::LocCtrl,"u0",u0);
00242 RegisterAt(base::LocCtrl,"TH0",TH0);
00243 RegisterAt(base::LocCtrl,"deltaT0",deltaT0);
00244 RegisterAt(base::LocCtrl,"L0",L0);
00245 RegisterAt(base::LocCtrl,"Re0",Re0);
00246 RegisterAt(base::LocCtrl,"Pr0",Pr0);
00247 RegisterAt(base::LocCtrl,"OutputErrorEvery",IntegrateEvery);
00248 }
00249 virtual void register_at(ControlDevice& Ctrl) {
00250 base::register_at(Ctrl);
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
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379 protected:
00380 double u0, TH0, deltaT0, L0, Re0, Pr0;
00381 interpolation_type* _Interpolation;
00382 int IntegrateEvery;
00383
00384 };