00001
00002
00003 #ifndef AMROC_PROBLEM_H
00004 #define AMROC_PROBLEM_H
00005
00006 #define DIM 3
00007
00008 #include "LBMProblem.h"
00009 #include "LBMD3Q19.h"
00010
00011 typedef double DataType;
00012 typedef LBMD3Q19<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 int direction=1;
00021
00022 class LBMSpecific : public LBMType {
00023 typedef LBMType base;
00024 public:
00025 LBMSpecific() : base(), omega(1.), force(0.), uf(0.2) {}
00026
00027 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00028 base::register_at(Ctrl,prefix);
00029 RegisterAt(base::LocCtrl,"Omega",omega);
00030 RegisterAt(base::LocCtrl,"uf",uf);
00031 RegisterAt(base::LocCtrl,"rho0",rho0);
00032 RegisterAt(base::LocCtrl,"u0",u0);
00033 RegisterAt(base::LocCtrl,"Direction",direction);
00034 }
00035
00036 virtual void SetupData(GridHierarchy* gh, const int& ghosts) {
00037 base::SetupData(gh,ghosts);
00038 if (base::LengthScale()!=DataType(1.) || base::TimeScale()!=DataType(1.)) {
00039 std::cerr << "dx(0) and dt(0) need to be equal to 1.0!" << std::endl;
00040 std::exit(-1);
00041 }
00042 if (omega<0.5 || omega>2.) {
00043 std::cerr << "0.5 <= omega <= 2 required." << std::endl;
00044 std::exit(-1);
00045 }
00046 base::SetGas(1.,base::LatticeViscosity(omega),base::LatticeSpeedOfSound());
00047 std::cout << "D3Q19: Gas_rho=" << base::GasDensity() << " Gas_Cs=" << base::GasSpeedofSound()
00048 << " Gas_nu=" << base::GasViscosity() << std::endl;
00049
00050 const double* bndry = base::GH().wholebndry();
00051 const int nbndry = base::GH().nbndry();
00052 DataType l0=bndry[nbndry*(1+2*1)]-bndry[nbndry*(0+2*1)];
00053 DataType uf=0.2;
00054 force=18.*base::LatticeViscosity(omega)*uf/l0/l0/l0;
00055 force = rho0*force/10.;
00056 if (direction<1) direction=1;
00057 if (direction>3) direction=3;
00058
00059 }
00060
00061 virtual double Step(vec_grid_data_type &fvec, vec_grid_data_type &ovec, vec_grid_data_type* Flux[],
00062 const double& t, const double& dt, const int& mpass) const {
00063 double ret = base::Step(fvec,ovec,Flux,t,dt,mpass);
00064
00065 int Nx = fvec.extents(0), Ny = fvec.extents(1), Nz = fvec.extents(2);
00066 int mxs = base::NGhosts(), mys = base::NGhosts(), mzs = base::NGhosts();
00067 int mxe = Nx-base::NGhosts()-1, mye = Ny-base::NGhosts()-1, mze = Nz-base::NGhosts()-1;
00068 base::MicroType *f = (base::MicroType *)fvec.databuffer();
00069
00070 for (register int k=mzs; k<=mze; k++)
00071 for (register int j=mys; j<=mye; j++)
00072 for (register int i=mxs; i<=mxe; i++) {
00073 if (direction==1) {
00074 f[base::idx(i,j,k,Nx,Ny)](1) += force;
00075 f[base::idx(i,j,k,Nx,Ny)](2) -= force;
00076 f[base::idx(i,j,k,Nx,Ny)](7) += force;
00077 f[base::idx(i,j,k,Nx,Ny)](8) -= force;
00078 f[base::idx(i,j,k,Nx,Ny)](9) += force;
00079 f[base::idx(i,j,k,Nx,Ny)](10) -= force;
00080 f[base::idx(i,j,k,Nx,Ny)](15) += force;
00081 f[base::idx(i,j,k,Nx,Ny)](16) -= force;
00082 f[base::idx(i,j,k,Nx,Ny)](17) += force;
00083 f[base::idx(i,j,k,Nx,Ny)](18) -= force;
00084 }
00085 else if (direction==2) {
00086 f[base::idx(i,j,k,Nx,Ny)](3) += force;
00087 f[base::idx(i,j,k,Nx,Ny)](4) -= force;
00088 f[base::idx(i,j,k,Nx,Ny)](7) += force;
00089 f[base::idx(i,j,k,Nx,Ny)](8) -= force;
00090 f[base::idx(i,j,k,Nx,Ny)](10) += force;
00091 f[base::idx(i,j,k,Nx,Ny)](9) -= force;
00092 f[base::idx(i,j,k,Nx,Ny)](11) += force;
00093 f[base::idx(i,j,k,Nx,Ny)](12) -= force;
00094 f[base::idx(i,j,k,Nx,Ny)](13) += force;
00095 f[base::idx(i,j,k,Nx,Ny)](14) -= force;
00096 }
00097 else if (direction==3) {
00098 f[base::idx(i,j,k,Nx,Ny)](5) += force;
00099 f[base::idx(i,j,k,Nx,Ny)](6) -= force;
00100 f[base::idx(i,j,k,Nx,Ny)](11) += force;
00101 f[base::idx(i,j,k,Nx,Ny)](12) -= force;
00102 f[base::idx(i,j,k,Nx,Ny)](14) += force;
00103 f[base::idx(i,j,k,Nx,Ny)](13) -= force;
00104 f[base::idx(i,j,k,Nx,Ny)](15) += force;
00105 f[base::idx(i,j,k,Nx,Ny)](16) -= force;
00106 f[base::idx(i,j,k,Nx,Ny)](18) += force;
00107 f[base::idx(i,j,k,Nx,Ny)](17) -= force;
00108 }
00109 }
00110
00111 return ret;
00112 }
00113
00114 private:
00115 DataType omega, force, uf;
00116 };
00117
00118
00119 class InitialConditionSpecific : public LBMInitialCondition<LBMType,DIM> {
00120 typedef LBMInitialCondition<LBMType,DIM> base;
00121 public:
00122 InitialConditionSpecific(LBMType &lbm) : base(lbm) {}
00123
00124 virtual void SetGrid(vec_grid_data_type& fvec, grid_data_type& workvec, const int& level) {
00125 DataType q[4] = { rho0, 0., 0., 0. };
00126 q[direction] = u0;
00127 LBM().ICStandard(fvec,LBMType::ConstantMacro,q,4,LBMType::Lattice);
00128 }
00129 };
00130
00131
00132 #endif