00001
00002
00003 #ifndef AMROC_PROBLEM_H
00004 #define AMROC_PROBLEM_H
00005
00006 #define DIM 2
00007
00008 #include "LBMProblem.h"
00009 #include "LBMD2Q9.h"
00010
00011 typedef double DataType;
00012 typedef LBMD2Q9<DataType> LBMType;
00013
00014 #define OWN_LBMSCHEME
00015 #define OWN_INITIALCONDITION
00016 #include "LBMStdProblem.h"
00017
00018 DataType rho0=1.;
00019 DataType u0=0.1;
00020
00021 class LBMSpecific : public LBMType {
00022 typedef LBMType base;
00023 public:
00024 LBMSpecific() : base(), omega(1.), force(0.), uf(0.2) {}
00025
00026 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00027 base::register_at(Ctrl,prefix);
00028 RegisterAt(base::LocCtrl,"Omega",omega);
00029 RegisterAt(base::LocCtrl,"uf",uf);
00030 RegisterAt(base::LocCtrl,"rho0",rho0);
00031 RegisterAt(base::LocCtrl,"u0",u0);
00032 }
00033
00034 virtual void SetupData(GridHierarchy* gh, const int& ghosts) {
00035 base::SetupData(gh,ghosts);
00036 if (base::LengthScale()!=DataType(1.) || base::TimeScale()!=DataType(1.)) {
00037 std::cerr << "dx(0) and dt(0) need to be equal to 1.0!" << std::endl;
00038 std::exit(-1);
00039 }
00040 if (omega<0.5 || omega>2.) {
00041 std::cerr << "0.5 <= omega <= 2 required." << std::endl;
00042 std::exit(-1);
00043 }
00044 base::SetGas(1.,base::LatticeViscosity(omega),base::LatticeSpeedOfSound());
00045 std::cout << "D2Q9: Gas_rho=" << base::GasDensity() << " Gas_Cs=" << base::GasSpeedofSound()
00046 << " Gas_nu=" << base::GasViscosity() << std::endl;
00047
00048 const double* bndry = base::GH().wholebndry();
00049 const int nbndry = base::GH().nbndry();
00050 DataType l0=bndry[nbndry*(1+2*1)]-bndry[nbndry*(0+2*1)];
00051 DataType uf=0.2;
00052 force=8.*base::LatticeViscosity(omega)*uf/l0/l0;
00053 force = rho0*force/6.;
00054 }
00055
00056 virtual double Step(vec_grid_data_type &fvec, vec_grid_data_type &ovec, vec_grid_data_type* Flux[],
00057 const double& t, const double& dt, const int& mpass) const {
00058 double ret = base::Step(fvec,ovec,Flux,t,dt,mpass);
00059
00060 int Nx = fvec.extents(0), Ny = fvec.extents(1);
00061 int mxs = base::NGhosts(), mys = base::NGhosts();
00062 int mxe = Nx-base::NGhosts()-1, mye = Ny-base::NGhosts()-1;
00063 base::MicroType *f = (base::MicroType *)fvec.databuffer();
00064
00065 for (register int j=mys; j<=mye; j++)
00066 for (register int i=mxs; i<=mxe; i++) {
00067 f[base::idx(i,j,Nx)](1) += force;
00068 f[base::idx(i,j,Nx)](4) += force;
00069 f[base::idx(i,j,Nx)](7) += force;
00070 f[base::idx(i,j,Nx)](2) -= force;
00071 f[base::idx(i,j,Nx)](5) -= force;
00072 f[base::idx(i,j,Nx)](8) -= force;
00073 }
00074
00075 return ret;
00076 }
00077
00078 private:
00079 DataType omega, force, uf;
00080 };
00081
00082
00083 class InitialConditionSpecific : public LBMInitialCondition<LBMType,DIM> {
00084 typedef LBMInitialCondition<LBMType,DIM> base;
00085 public:
00086 InitialConditionSpecific(LBMType &lbm) : base(lbm) {}
00087
00088 virtual void SetGrid(vec_grid_data_type& fvec, grid_data_type& workvec, const int& level) {
00089 DataType q[3] = { rho0, u0, 0. };
00090 LBM().ICStandard(fvec,LBMType::ConstantMacro,q,3,LBMType::Lattice);
00091 }
00092 };
00093
00094
00095 #endif