00001
00002
00003 #ifndef AMROC_PROBLEM_H
00004 #define AMROC_PROBLEM_H
00005
00006 #define DIM 1
00007
00008 #include "LBMProblem.h"
00009 #include "LBMD1Q3.h"
00010
00011 typedef double DataType;
00012 typedef LBMD1Q3<DataType> LBMType;
00013
00014 #define OWN_LBMSCHEME
00015 #define OWN_INITIALCONDITION
00016 #include "LBMStdProblem.h"
00017
00018 DataType rho0=0.5;
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.0) {}
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 << "D1Q3: 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.0;
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);
00061 int mxs = base::NGhosts();
00062 int mxe = Nx-base::NGhosts()-1;
00063 base::MicroType *f = (base::MicroType *)fvec.databuffer();
00064
00065 for (register int i=mxs; i<=mxe; i++) {
00066 f[i](1) += force;
00067 f[i](2) -= force;
00068 }
00069
00070 return ret;
00071 }
00072
00073 private:
00074 DataType omega, force, uf;
00075 };
00076
00077
00078 class InitialConditionSpecific : public LBMInitialCondition<LBMType,DIM> {
00079 typedef LBMInitialCondition<LBMType,DIM> base;
00080 public:
00081 InitialConditionSpecific(LBMType &lbm) : base(lbm) {}
00082
00083 DataType caux[6];
00084
00085 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00086 base::LocCtrl = Ctrl.getSubDevice(prefix+"InitialConditionCustom");
00087 RegisterAt(base::LocCtrl,"Type",Type);
00088 RegisterAt(base::LocCtrl,"Scaling",Scaling);
00089 RegisterAt(base::LocCtrl,"NAux",NAux);
00090 char VariableName[16];
00091
00092 for (int i=0; i<LBM().NMicroVar(); i++) {
00093 caux[i]=DataType(0.);
00094 std::sprintf(VariableName,"CAux(%d)",i+1);
00095 RegisterAt(base::LocCtrl,VariableName,caux[i]);
00096 }
00097 for (int i=LBM().NMicroVar(); i<6; i++) {
00098 caux[i]=DataType(1.);
00099 std::sprintf(VariableName,"CAux(%d)",i+1);
00100 RegisterAt(base::LocCtrl,VariableName,caux[i]);
00101 }
00102 for (int i=0; i<5; i++) {
00103 std::cout << "CAux(" << i << ")=" << caux[i] <<std::endl;
00104 }
00105
00106 }
00107
00108 virtual void SetGrid(vec_grid_data_type& fvec, grid_data_type& workvec, const int& level) {
00109 BBox box = fvec.bbox();
00110 int mxs = box.lower(0)/box.stepsize(0), mxe = box.upper(0)/box.stepsize(0);
00111 int b = fvec.bottom();
00112
00113 for (int i=0; i<5; i++) {
00114 std::cout << "SetGrid CAux(" << i << ")=" << caux[i] <<std::endl;
00115 }
00116 DataType a = caux[2];
00117 DataType c = caux[3];
00118 DataType w = caux[4]/ DataType(2.);
00119 base::MicroType *f = (base::MicroType *)fvec.databuffer();
00120 base::MacroType q;
00121 q(1) = caux[1];
00122 for (register int i=mxs; i<=mxe; i++) {
00123 if (fabs(i-c)<=w) q(0) = a*caux[0];
00124 else q(0) = caux[0];
00125 f[b+i] = LBM().Equilibrium(q);
00126 }
00127 }
00128 };
00129
00130
00131 #endif