fsi/motion-amroc/HingedWing/src/FluidProblem.h

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 
00003 #ifndef AMROC_FLUIDPROBLEM_H
00004 #define AMROC_FLUIDPROBLEM_H
00005 #define DIM 3
00006 
00007 #include "LBMProblem.h"
00008 #include "LBMD3Q19.h"
00009 
00010 typedef double DataType;
00011 typedef LBMD3Q19<DataType> LBMType;
00012 
00013 #define OWN_LBMSCHEME
00014 #define OWN_AMRSOLVER
00015 //#define OWN_INITIALCONDITION
00016 //#define OWN_BOUNDARYCONDITION
00017 #include "LBMStdProblem.h"
00018 #include "F77Interfaces/F77GFMBoundary.h"
00019 #include "AMRELCGFMSolver.h"
00020 #include "LBMGFMBoundary.h"
00021 #include "Interfaces/SchemeGFMFileOutput.h"
00022 
00023 class LBMSpecific : public LBMType {
00024   typedef LBMType base;
00025   typedef ads::FixedArray<8,DataType> PVType;
00026   typedef ads::FixedArray<10,DataType> PVCType;
00027 public:
00028   LBMSpecific() : base(), omega(1.), u_p_avg(1.), l0_p(-1.) {
00029 
00030   }
00031 
00032   virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00033     base::register_at(Ctrl,prefix);
00034     RegisterAt(base::LocCtrl,"Omega",omega);
00035     RegisterAt(base::LocCtrl,"u_p_avg",u_p_avg);
00036     RegisterAt(base::LocCtrl,"l0_p",l0_p);
00037 
00038   }
00039 
00040   virtual void SetupData(GridHierarchy* gh, const int& ghosts) {
00041     base::SetupData(gh,ghosts);
00042 
00043     const double* bndry = base::GH().wholebndry();
00044     const int nbndry = base::GH().nbndry();
00045     if (l0_p <= 0. ) {
00046       l0_p=bndry[nbndry*(1+2*1)]-bndry[nbndry*(0+2*1)];
00047       if (nbndry>1) l0_p=bndry[nbndry*(1+2*1)]-bndry[1+nbndry*(1+2*1)];
00048 
00049     }
00050     DataType omega = base::Omega(base::TimeScale());
00051     //DataType u_p_avg=Side(0).aux[0];
00052     //DataType u_p_avg=WIND;//caux[1];
00053     double Re_p=u_p_avg*l0_p/base::GasViscosity();
00054     double Re_lbm=(u_p_avg/base::VelocityScale())*(l0_p/base::LengthScale())/base::LatticeViscosity(omega);
00055     ( comm_service::log() << "l0_p=" << l0_p << "   u_p_avg=" << u_p_avg
00056                           << "   Re_p=" << Re_p << "   Re_lbm=" << Re_lbm
00057                           << "\n   omega=" << omega
00058                           << "   SpeedUp=" << base::SpeedUp()
00059                           << "   u_avg_LBM=" << u_p_avg/base::VelocityScale()
00060                           << "   SpeedRatio=" << u_p_avg/base::VelocityScale()*std::sqrt(DataType(3.))
00061                           << "   base::TO()=" << base::T0()
00062                           << "\n   Re_p-Re_lbm=" << Re_p-Re_lbm
00063                           << std::endl ).flush();
00064     //assert(std::fabs(Re_p-Re_lbm)<DBL_EPSILON*LB_FACTOR*LB_FACTOR*LB_FACTOR);
00065     if ( std::fabs(Re_p-Re_lbm)>DBL_EPSILON*LB_FACTOR*LB_FACTOR*LB_FACTOR ) {
00066       std::cerr << "Physical scaling requires fabs(Re_p-Re_lbm)<DBL_EPSILON*LB_FACTOR*LB_FACTOR*LB_FACTOR !" << std::endl;
00067       std::cerr << "Physical scaling requires fabs(Re_p-Re_lbm)<DBL_EPSILON*LB_FACTOR*LB_FACTOR*LB_FACTOR !" << std::endl;
00068       std::cerr << "Re_p " << Re_p << " Re_lbm " << Re_lbm << " std::fabs(Re_p-Re_lbm) " << std::fabs(Re_p-Re_lbm) << std::endl;
00069       std::exit(-1);
00070     }
00071     if (omega<0.5 || omega>2.) {
00072       std::cerr << "0.5 <= omega <= 2 required." << std::endl;
00073       std::exit(-1);
00074     }
00075     base::SetGas(base::GasDensity(),base::GasViscosity(),base::GasSpeedofSound());
00076     ( comm_service::log() << "FluidProblem SetupData() D3Q19: Gas_rho=" << base::GasDensity() << "   Gas_Cs=" << base::GasSpeedofSound()
00077                           << "   Gas_nu=" << base::GasViscosity() << std::endl ).flush();
00078 
00079   }
00080 
00081   virtual double Step(vec_grid_data_type &fvec, vec_grid_data_type &ovec, vec_grid_data_type* Flux[],
00082                       const double& t, const double& dt, const int& mpass) const {
00083 
00084     LBMType::Step(fvec, ovec, Flux, t, dt, mpass);
00085     return DataType(1.);
00086   }
00087 
00088 private:
00089   DataType omega, u_p_avg, l0_p;
00090 
00091 protected:
00092 
00093 };
00094 
00095 /*class InitialConditionSpecific : public LBMInitialCondition<LBMType,DIM> {
00096   typedef LBMInitialCondition<LBMType,DIM> base;
00097 public:
00098   InitialConditionSpecific(LBMType &lbm) : base(lbm) {}
00099   typedef base::MicroType MicroType;
00100   typedef base::MacroType MacroType;
00101 
00102   virtual void SetGrid(vec_grid_data_type& fvec, grid_data_type& workvec, const int& level) {
00103     MacroType q, q_;
00104     DataType elevation = 15., profile, height, ldx;
00105     BBox wb = base::GH().wholebbox();
00106     BBox bb = fvec.bbox();
00107     DCoords wc;
00108 
00109     q_(0) = aux[0]/LBM().DensityScale();
00110     q_(1) = aux[1]/LBM().VelocityScale();
00111     q_(2) = aux[2]/LBM().VelocityScale();
00112     q_(3) = aux[3]/LBM().VelocityScale();
00113 
00114     MicroType *f = (MicroType *)fvec.databuffer();
00115 
00116     int Nx = fvec.extents(0), Ny = fvec.extents(1);
00117     int b = fvec.bottom();
00118     assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1));
00119     int mxs = std::max(bb.lower(0)/bb.stepsize(0),fvec.lower(0)/fvec.stepsize(0));
00120     int mys = std::max(bb.lower(1)/bb.stepsize(1),fvec.lower(1)/fvec.stepsize(1));
00121     int mzs = std::max(bb.lower(2)/bb.stepsize(2),fvec.lower(2)/fvec.stepsize(2));
00122     int mxe = std::min(bb.upper(0)/bb.stepsize(0),fvec.upper(0)/fvec.stepsize(0));
00123     int mye = std::min(bb.upper(1)/bb.stepsize(1),fvec.upper(1)/fvec.stepsize(1));
00124     int mze = std::min(bb.upper(2)/bb.stepsize(2),fvec.upper(2)/fvec.stepsize(2));
00125     int lc_indx[3];
00126 
00127     lc_indx[0] = base::GH().worldCoords(fvec.lower(),fvec.stepsize())(0);
00128     lc_indx[1] = base::GH().worldCoords(fvec.lower(),fvec.stepsize())(1);
00129 
00130     ldx = DataType(fvec.stepsize(2))/DataType(wb.stepsize(2))*LBM().L0();
00131 
00132     for (register int k=mzs; k<=mze; k++) {
00133       lc_indx[2] = (int) std::floor(DataType(k)*DataType(ldx));
00134       wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
00135       height = DataType(wc(2))*DataType(wb.stepsize(2));
00136 
00137       //std::cout << " height " << height << std::endl;
00138       if ( height < elevation ) {
00139         profile = DataType(height*height)/DataType(elevation*elevation);
00140 
00141         q(0) = q_(0);
00142         q(1) = q_(1)*profile;
00143         q(2) = q_(2)*profile;
00144         q(3) = q_(3)*profile;
00145       }
00146       else {
00147         q = q_;
00148 
00149       }
00150       for (register int j=mys; j<=mye; j++)
00151         for (register int i=mxs; i<=mxe; i++)
00152           f[b+LBM().idx(i,j,k,Nx,Ny)] = LBM().Equilibrium(q);
00153     }
00154 
00155   }
00156 };
00157 
00158 class BoundaryConditionsSpecific : public LBMBoundaryConditions<LBMType,DIM> {
00159   typedef LBMBoundaryConditions<LBMType,DIM> base;
00160 public:
00161   typedef base::MicroType MicroType;
00162   typedef base::MacroType MacroType;
00163   BoundaryConditionsSpecific(LBMType &lbm) : base(lbm) {}
00164 
00165   DataType caux[6];
00166   int ctype, cside, cscaling, ncaux;
00167 
00168   virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00169     base::register_at(Ctrl, prefix);
00170     base::LocCtrl = Ctrl.getSubDevice(prefix+"BoundaryConditionCustom");
00171     RegisterAt(base::LocCtrl,"Type",ctype);
00172     RegisterAt(base::LocCtrl,"Side",cside);
00173     RegisterAt(base::LocCtrl,"Scaling",cscaling);
00174     RegisterAt(base::LocCtrl,"NCAux",ncaux);
00175     char VariableName[16];
00176     for (int i=0; i<6; i++) {
00177       caux[i]=DataType(1.);
00178       std::sprintf(VariableName,"CAux(%d)",i+1);
00179       RegisterAt(base::LocCtrl,VariableName,caux[i]);
00180     }
00181   }
00182 
00183   virtual void SetupData(GridHierarchy* gh, const int& ghosts) {
00184     base::SetupData(gh,ghosts);
00185     const double* bndry = base::GH().wholebndry();
00186     const int nbndry = base::GH().nbndry();
00187     DataType l0_p=bndry[nbndry*(1+2*1)]-bndry[nbndry*(0+2*1)];
00188     if (nbndry>1) l0_p=bndry[nbndry*(1+2*1)]-bndry[1+nbndry*(1+2*1)];
00189     DataType omega = LBM().Omega(LBM().TimeScale());
00190     //DataType u_p_avg=Side(0).aux[0];
00191     DataType u_p_avg=caux[1];
00192     double Re_p=u_p_avg*l0_p/LBM().GasViscosity();
00193     double Re_lbm=(u_p_avg/LBM().VelocityScale())*(l0_p/LBM().LengthScale())/LBM().LatticeViscosity(omega);
00194     ( comm_service::log() << "l0_p=" << l0_p << "   u_p_avg=" << u_p_avg
00195               << "   Re_p=" << Re_p << "   Re_lbm=" << Re_lbm
00196               << "\n   omega=" << omega
00197               << "   SpeedUp=" << LBM().SpeedUp()
00198               << "   u_avg_LBM=" << u_p_avg/LBM().VelocityScale()
00199               << "   SpeedRatio=" << u_p_avg/LBM().VelocityScale()*std::sqrt(DataType(3.))
00200               << "   LBM().TO()=" << LBM().T0()
00201               << "\n   Re_p-Re_lbm=" << Re_p-Re_lbm
00202               << std::endl ).flush();
00203     assert(std::fabs(Re_p-Re_lbm)<DBL_EPSILON*LB_FACTOR*LB_FACTOR*LB_FACTOR);
00204   }
00205 
00206   virtual void SetBndry(vec_grid_data_type &fvec, const int& level, const BBox &bb,
00207                         const int &dir, const double& time) {
00208 
00209     //std::cout << " fluid problem custom boundary begin\n";
00210 
00211     base::SetBndry(fvec,level,bb,dir,time);
00212     BBox wholebbox = base::GH().wholebbox();
00213     BBox inside = bb*coarsen(wholebbox,wholebbox.stepsize()/wholebbox.stepsize());
00214     if (!inside.empty())
00215       LBM().BCStandard(fvec,bb,LBMType::NoSlipWall,dir);
00216 
00217     int scaling = LBMType::Physical;
00218     int Nx = fvec.extents(0), Ny = fvec.extents(1);//, Nz = fvec.extents(2);
00219     int b = fvec.bottom();
00220     assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1));
00221     int mys = std::max(bb.lower(1)/bb.stepsize(1),fvec.lower(1)/fvec.stepsize(1));
00222     int mzs = std::max(bb.lower(2)/bb.stepsize(2),fvec.lower(2)/fvec.stepsize(2));
00223     int mxe = std::min(bb.upper(0)/bb.stepsize(0),fvec.upper(0)/fvec.stepsize(0));
00224     int mye = std::min(bb.upper(1)/bb.stepsize(1),fvec.upper(1)/fvec.stepsize(1));
00225     int mze = std::min(bb.upper(2)/bb.stepsize(2),fvec.upper(2)/fvec.stepsize(2));
00226     MicroType *f = (MicroType *)fvec.databuffer();
00227 
00228     DataType rho=caux[0], a0=caux[1], a1=caux[2], a2=caux[3], elevation=caux[4], jr, hr;
00229 
00230     MacroType q, q_;
00231     DataType profile, height, ldx;
00232 
00233     if (scaling==LBMType::Physical) {
00234       rho /= LBM().DensityScale();
00235       a0 /= LBM().VelocityScale();
00236       a1 /= LBM().VelocityScale();
00237       a2 /= LBM().VelocityScale();
00238       elevation /= LBM().L0();
00239     }
00240 
00241     if (bb.upper(0)+bb.stepsize(0)==wholebbox.lower(0)) {
00242 
00243 
00244       BBox wb = base::GH().wholebbox();
00245       DCoords wc;
00246       q_(0) = rho;
00247       q_(1) = a0;
00248       q_(2) = a1;///LBM().VelocityScale(); // LBM().SpeedUp()*aux[2]/LBM().VelocityScale();
00249       q_(3) = a2;///LBM().VelocityScale(); // LBM().SpeedUp()*aux[3]/LBM().VelocityScale();
00250 
00251       int lc_indx[3];
00252 
00253       lc_indx[0] = base::GH().worldCoords(fvec.lower(),fvec.stepsize())(0);
00254       lc_indx[1] = base::GH().worldCoords(fvec.lower(),fvec.stepsize())(1);
00255       ldx = DataType(fvec.stepsize(2))/DataType(wb.stepsize(2))*DataType(LBM().L0());
00256 
00257       for (register int k=mzs; k<=mze; k++) {
00258         lc_indx[2] = k;
00259 
00260         wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
00261         height = DataType(wc(2))*DataType(wb.stepsize(2));
00262         if ( height < elevation ) {
00263           profile = DataType(height*height)/DataType(elevation*elevation);
00264 
00265           q(0) = q_(0);
00266           q(1) = q_(1)*profile;//(k*LBM().L0())*(k*LBM().L0())/(elevation*elevation);
00267           q(2) = q_(2)*profile;//(k*LBM().L0())/elevation;
00268           q(3) = q_(3)*profile;//(k*LBM().L0())/elevation;
00269         }
00270         else {
00271           q = q_;
00272           //std::cout << " level " << level << " mzs " << mzs <<" k " << k << " mze " << mze << "\t height " << height << " elevation " << elevation << " q " << q << " ldx " << ldx << std::endl;
00273 
00274         }
00275         for (register int j=mys; j<=mye; j++) {
00276           f[b+LBM().idx(mxe,j,k,Nx,Ny)] = LBM().Equilibrium(q);
00277         }
00278 
00279 
00280       }
00281     }
00282   }
00283 };*/
00284 
00285 
00286 template <class LBMType, int dim>
00287 class LBMELCGFMBoundary : public LBMGFMBoundary<LBMType,dim> {
00288   typedef typename LBMType::MicroType MicroType;
00289   typedef typename MicroType::InternalDataType DataType;
00290   typedef LBMGFMBoundary<LBMType,dim> base;
00291 public:
00292   typedef typename base::vec_grid_data_type vec_grid_data_type;
00293   typedef typename base::grid_data_type grid_data_type;
00294   typedef typename base::point_type point_type;
00295   typedef AMRELCGFMSolver<MicroType,FixupType,FlagType,dim> amr_elc_solver;
00296 
00297   LBMELCGFMBoundary(LBMType &scheme, amr_elc_solver& solver) : base(scheme), _solver(solver) {
00298     base::SetNAux(base::Dim());
00299   }
00300 
00301   virtual void SetBndryAux(vec_grid_data_type& gdu, grid_data_type& gdphi,
00302                            const MicroType* u, DataType* aux, const int& Level,
00303                            double t, const int& nc, const int* idx,
00304                            const point_type* xc, const DataType* distance,
00305                            const point_type* normal) {
00306     _solver.ExtractBoundaryVelocities(gdu.bbox(),nc,xc,reinterpret_cast<point_type *>(aux));
00307   }
00308 
00309 protected:
00310   amr_elc_solver& _solver;
00311 };
00312 
00313 class SolverObjects {
00314 public:
00315   SolverObjects() {
00316     _LBMSpecific = new LBMSpecific();
00317     _IntegratorSpecific = new IntegratorSpecific(*_LBMSpecific);
00318     _InitialConditionSpecific = new InitialConditionSpecific(*_LBMSpecific);
00319     _BoundaryConditionsSpecific = new BoundaryConditionsSpecific(*_LBMSpecific);
00320   }
00321 
00322   ~SolverObjects() {
00323     delete _BoundaryConditionsSpecific;
00324     delete _InitialConditionSpecific;
00325     delete _IntegratorSpecific;
00326     delete _LBMSpecific;
00327   }
00328 
00329 protected:
00330   LBMSpecific *_LBMSpecific;
00331   IntegratorSpecific *_IntegratorSpecific;
00332   InitialConditionSpecific *_InitialConditionSpecific;
00333   BoundaryConditionsSpecific *_BoundaryConditionsSpecific;
00334 };
00335 
00336 class FluidSolverSpecific : protected SolverObjects,
00337     public AMRELCGFMSolver<MicroType,FixupType,FlagType,DIM> {
00338   typedef MicroType::InternalDataType DataType;
00339   typedef AMRELCGFMSolver<MicroType,FixupType,FlagType,DIM> base;
00340   typedef AMRInterpolation<MicroType,DIM> interpolation_type;
00341   typedef interpolation_type::point_type point_type;
00342   typedef SchemeGFMFileOutput<LBMType,FixupType,FlagType,DIM> output_type;
00343   typedef ads::FixedArray<10,DataType> PVCType;
00344 public:
00345   FluidSolverSpecific() : SolverObjects(), base(*SolverObjects::_IntegratorSpecific,
00346                                                 *SolverObjects::_InitialConditionSpecific,
00347                                                 *SolverObjects::_BoundaryConditionsSpecific) {
00348     SetLevelTransfer(new LBMF77LevelTransfer<LBMType,DIM>(SolverObjects::_IntegratorSpecific->Scheme(), f_prolong, f_restrict));
00349     SetFileOutput(new output_type(*this,SolverObjects::_IntegratorSpecific->Scheme()));
00350 
00351     SetFixup(new FixupSpecific(SolverObjects::_IntegratorSpecific->Scheme()));
00352     SetFlagging(new FlaggingSpecific(*this,SolverObjects::_IntegratorSpecific->Scheme()));
00353     AddGFM(new GhostFluidMethod<MicroType,DIM>(new LBMELCGFMBoundary<LBMType,DIM>(SolverObjects::_IntegratorSpecific->LBM(),*this),
00354                                                new CPTLevelSet<DataType,DIM>()));
00355     //ABL_CPT = new CPTLevelSet<DataType,DIM>();
00356     SetCoupleGFM(0);
00357 
00358   }
00359 
00360   ~FluidSolverSpecific() {
00361     DeleteGFM(_GFM[0]);
00362     delete _LevelTransfer;
00363     delete _Flagging;
00364     delete _Fixup;
00365     delete _FileOutput;
00366     delete _BoundaryConditionsSpecific;
00367     _BoundaryConditionsSpecific = (BoundaryConditionsSpecific *) 0;
00368     delete _InitialConditionSpecific;
00369     _InitialConditionSpecific = (InitialConditionSpecific *) 0;
00370     delete _IntegratorSpecific;
00371     _IntegratorSpecific = (IntegratorSpecific *) 0;
00372     delete _Interpolation;
00373     _Interpolation = (AMRELCGFMSolver<MicroType,FixupType,FlagType,DIM>::interpolation_type *) 0;
00374   }
00375 
00376   virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00377     base::register_at(Ctrl,prefix);
00378 
00379   }
00380   virtual void register_at(ControlDevice& Ctrl) {
00381     base::register_at(Ctrl);
00382   }
00383 
00384   virtual void SendBoundaryData() {
00385     START_WATCH
00386         for (register int l=0; l<=FineLevel(base::GH()); l++) {
00387       int Time = CurrentTime(base::GH(),l);
00388       ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time,
00389                                               l, base::Dim()+4, base::t[l]);
00390       if (CurrentTime(base::GH(),base::CouplingLevel) != Time)
00391         ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time+TimeStep(base::U(),l),
00392                                                 l, base::Dim()+4, base::t[l]+base::dt[l]);
00393     }
00394     END_WATCH(FLUID_CPL_PRESSURE_CALCULATE)
00395         base::SendBoundaryData();
00396   }
00397 
00398   virtual void SetupData() {
00399     ( comm_service::log() << " fluid setup begun" << std::endl ).flush();
00400     base::SetupData();
00401     base::AdaptBndTimeInterpolate = 0;
00402     base::NAMRTimeSteps = 1;
00403     base::Step[0].LastTime = 1.e37;
00404     base::Step[0].VariableTimeStepping = -1;
00405     base::Step[0].dtv[0] = ((LBMIntegrator<LBMType,DIM> &)Integrator_()).LBM().TimeScale();
00406 
00407     _Interpolation->SetupData(base::PGH(), base::NGhosts());
00408 
00409   }
00410 
00411 protected:
00412 
00413 };