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_AMRSOLVER
00015 #include "LBMStdProblem.h"
00016
00017 class SolverSpecific : public AMRSolver<MicroType,FixupType,FlagType,DIM> {
00018 typedef AMRSolver<MicroType,FixupType,FlagType,DIM> base;
00019 typedef MicroType::InternalDataType DataType;
00020 public:
00021 SolverSpecific(IntegratorSpecific& integ,
00022 base::initial_condition_type& init,
00023 base::boundary_conditions_type& bc) : base(integ, init, bc) {
00024 SetLevelTransfer(new LBMF77LevelTransfer<LBMType,DIM>(integ.Scheme(), f_prolong, f_restrict));
00025 SetFileOutput(new SchemeFileOutput<LBMType,DIM>(integ.Scheme()));
00026 SetFixup(new FixupSpecific(integ.Scheme()));
00027 SetFlagging(new FlaggingSpecific(*this,integ.Scheme()));
00028 }
00029
00030 ~SolverSpecific() {
00031 delete _LevelTransfer;
00032 delete _Flagging;
00033 delete _Fixup;
00034 delete _FileOutput;
00035 }
00036
00037 virtual void SetupData() {
00038 base::SetupData();
00039 base::NAMRTimeSteps = 1;
00040 base::AdaptBndTimeInterpolate = 0;
00041 base::Step[0].LastTime = 1.e37;
00042 base::Step[0].VariableTimeStepping = -1;
00043 base::Step[0].dtv[0] = ((LBMIntegrator<LBMType,DIM> &)Integrator_()).LBM().TimeScale();
00044
00045 LBMType &lbm = (LBMType &)((LBMIntegrator<LBMType,DIM> &)Integrator_()).LBM();
00046 DataType u_avg_p = ((LBMBoundaryConditions<LBMType,DIM> &)BoundaryConditions_()).Side(3).aux[0];
00047 BBox wb = base::GH().wholebbox();
00048 DCoords lc = base::GH().worldCoords(wb.lower(),wb.stepsize());
00049 DCoords uc = base::GH().worldCoords(wb.upper()+wb.stepsize(),wb.stepsize());
00050 DataType l0_p=uc(0)-lc(0);
00051 DataType omega = lbm.Omega(lbm.TimeScale());
00052 double Re_p=u_avg_p*l0_p/lbm.GasViscosity();
00053 double Re_lbm=(u_avg_p/lbm.VelocityScale())*(l0_p/lbm.LengthScale())/lbm.LatticeViscosity(omega);
00054 std::cout << "l0_p=" << l0_p << " u_avg_p=" << u_avg_p
00055 << " Re_p=" << Re_p << " Re_lbm=" << Re_lbm << std::endl;
00056 if ( lbm.TurbulenceType() == 1 ) {
00057 std::cout << "LES_Smagorinsky : Smagorinsky Constant = " << lbm.SmagorinskyConstant()
00058 << std::endl;
00059 }
00060 assert(std::fabs(Re_p-Re_lbm)<DBL_EPSILON*LB_FACTOR);
00061 }
00062 };