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 #define OWN_AMRSOLVER
00017 #include "LBMStdProblem.h"
00018 #include "AMRGFMSolver.h"
00019 #include "LBMGFMBoundary.h"
00020 #include "Interfaces/SchemeGFMFileOutput.h"
00021
00022 DataType rho0=1.;
00023 DataType rd, xs[2], rs, T, Tr;
00024 DataType pi=4.0*std::atan(1.0);
00025
00026 class LBMSpecific : public LBMType {
00027 typedef LBMType base;
00028 public:
00029 LBMSpecific() : base(), omega(1.) {}
00030
00031 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00032 base::register_at(Ctrl,prefix);
00033 RegisterAt(base::LocCtrl,"Omega",omega);
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 << "D2Q9: Gas_rho=" << base::GasDensity() << " Gas_Cs=" << base::GasSpeedofSound()
00048 << " Gas_nu=" << base::GasViscosity() << std::endl;
00049 }
00050
00051 private:
00052 DataType omega;
00053 };
00054
00055
00056 class InitialConditionSpecific : public LBMInitialCondition<LBMType,DIM> {
00057 typedef LBMInitialCondition<LBMType,DIM> base;
00058 public:
00059 InitialConditionSpecific(LBMType &lbm) : base(lbm) {}
00060
00061 virtual void SetGrid(vec_grid_data_type& fvec, grid_data_type& workvec, const int& level) {
00062 DataType q[3] = { rho0, 0., 0. };
00063 LBM().ICStandard(fvec,LBMType::ConstantMacro,q,3,LBMType::Lattice);
00064 }
00065 };
00066
00067
00068
00069 class GFMBoundarySpecific : public LBMGFMBoundary<LBMType,DIM> {
00070 typedef LBMGFMBoundary<LBMType,DIM> base;
00071 public:
00072 GFMBoundarySpecific(LBMType &lbm) : base(lbm) {}
00073
00074 virtual void SetBndry(vec_grid_data_type& gdu, const int& nc, const int* idx,
00075 const MicroType* u, const point_type* xc, const DataType* distance,
00076 const point_type* normal, DataType* aux, const int naux, const double t)
00077 { LBM().GFMBCStandard(gdu,base::Type,nc,idx,u,xc,distance,normal,aux,naux,LBMType::Lattice); }
00078
00079 virtual void SetBndryAux(vec_grid_data_type& gdu, grid_data_type& gdphi,
00080 const MicroType* u, DataType* aux, const int& Level,
00081 double t, const int& nc, const int* idx,
00082 const point_type* xc, const DataType* distance,
00083 const point_type* normal) {
00084 DataType xsc[2], vc[2];
00085 xsc[0] = xs[0]; xsc[1] = xs[1];
00086 vc[0] = 0.; vc[1] = 0.;
00087 if (T!=0.) {
00088 xsc[0] += rs*std::cos(2.*pi/T*t);
00089 xsc[1] += rs*std::sin(2.*pi/T*t);
00090 vc[0] = -2.*pi/T*std::sin(2.*pi/T*t);
00091 vc[1] = 2.*pi/T*std::cos(2.*pi/T*t);
00092 }
00093 for (register int i=0; i<nc; i++) {
00094
00095
00096 DataType alpha = std::atan((xc[i](1)-xsc[1])/(xc[i](0)-xsc[0]));
00097 if (xc[i](0)<xsc[0]) alpha += pi;
00098 DataType xp = xsc[0]+rd*std::cos(alpha)-xs[0];
00099 DataType yp = xsc[1]+rd*std::sin(alpha)-xs[1];
00100 DataType rc = std::sqrt(xp*xp+yp*yp);
00101 aux[i*base::NAux()] = rc*vc[0];
00102 aux[i*base::NAux()+1] = rc*vc[1];
00103 if (Tr!=0.) {
00104 DataType vr = 2.*pi/Tr*rd;
00105 aux[i*base::NAux()] -= vr*std::sin(alpha);
00106 aux[i*base::NAux()+1] += vr*std::cos(alpha);
00107 }
00108 }
00109 }
00110 };
00111
00112
00113 class GFMLevelSetSpecific : public GFMLevelSet<DataType,DIM> {
00114 typedef GFMLevelSet<DataType,DIM> base;
00115
00116 public:
00117 typedef base::grid_fct_type grid_fct_type;
00118 typedef base::grid_data_type grid_data_type;
00119
00120 GFMLevelSetSpecific() {}
00121 virtual ~GFMLevelSetSpecific() {}
00122
00123 virtual void SetGrid(grid_data_type& gdphi, const int& Level, const double& t) {
00124 int Nx = gdphi.extents(0), Ny = gdphi.extents(1);
00125 DCoords lb = base::GH().worldCoords(gdphi.lower(), gdphi.stepsize());
00126 DCoords dx = base::GH().worldStep(gdphi.stepsize());
00127 DataType *phi = (DataType *)gdphi.databuffer();
00128
00129 DataType xsc[2];
00130 xsc[0] = xs[0]; xsc[1] = xs[1];
00131 if (T!=0.) {
00132 xsc[0] += rs*std::cos(2.*pi/T*t);
00133 xsc[1] += rs*std::sin(2.*pi/T*t);
00134 }
00135 for (register int j=0; j<Ny; j++) {
00136 DataType yd = (DataType(j)+static_cast<DataType>(0.5))*dx(1) + lb(1) - xsc[1];
00137 for (register int i=0; i<Nx; i++) {
00138 DataType xd = (DataType(i)+static_cast<DataType>(0.5))*dx(0) + lb(0) - xsc[0];
00139 phi[i+Nx*j] = std::sqrt(xd*xd+yd*yd) - rd;
00140 }
00141 }
00142 }
00143 };
00144
00145
00146 class GhostFluidMethodSpecific : public GhostFluidMethod<MicroType,DIM> {
00147 typedef GhostFluidMethod<MicroType,DIM> base;
00148 public:
00149 GhostFluidMethodSpecific(boundary_type* boundary, levelset_type* levelset) :
00150 base(boundary,levelset) {
00151 xs[0]=0.; xs[1]=0.; rs = 1.; rd=1.; T=100.; Tr=0.;
00152 }
00153
00154 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00155 base::register_at(Ctrl,prefix);
00156 RegisterAt(base::LocCtrl,"Center(1)",xs[0]);
00157 RegisterAt(base::LocCtrl,"Center(2)",xs[1]);
00158 RegisterAt(base::LocCtrl,"RotationDuration",Tr);
00159 RegisterAt(base::LocCtrl,"CircleDuration",T);
00160 RegisterAt(base::LocCtrl,"RadiusCircle",rs);
00161 RegisterAt(base::LocCtrl,"RadiusSphere",rd);
00162 }
00163 };
00164
00165
00166 class SolverSpecific :
00167 public AMRGFMSolver<MicroType,FixupType,FlagType,DIM> {
00168 typedef AMRGFMSolver<MicroType,FixupType,FlagType,DIM> base;
00169 public:
00170 SolverSpecific(IntegratorSpecific& integ,
00171 base::initial_condition_type& init,
00172 base::boundary_conditions_type& bc) : base(integ, init, bc) {
00173 SetLevelTransfer(new LBMF77LevelTransfer<LBMType,DIM>(integ.Scheme(), f_prolong, f_restrict));
00174 SetFileOutput(new SchemeGFMFileOutput<LBMType,FixupType,FlagType,DIM>(*this,integ.Scheme()));
00175 SetFixup(new FixupSpecific(integ.Scheme()));
00176 SetFlagging(new FlaggingSpecific(*this,integ.Scheme()));
00177 AddGFM(new GhostFluidMethodSpecific(new GFMBoundarySpecific(integ.LBM()),
00178 new GFMLevelSetSpecific()));
00179 }
00180
00181 ~SolverSpecific() {
00182 DeleteGFM(_GFM[0]);
00183 delete _LevelTransfer;
00184 delete _Flagging;
00185 delete _Fixup;
00186 delete _FileOutput;
00187 }
00188
00189 virtual void SetupData() {
00190 base::SetupData();
00191 base::NAMRTimeSteps = 1;
00192 base::AdaptBndTimeInterpolate = 0;
00193 base::Step[0].LastTime = 1.e37;
00194 base::Step[0].VariableTimeStepping = -1;
00195 base::Step[0].dtv[0] = ((LBMIntegrator<LBMType,DIM> &)Integrator_()).LBM().TimeScale();
00196 }
00197 };
00198
00199 #endif