00001
00002
00003
00004 #ifndef AMROC_PROBLEM_H
00005 #define AMROC_PROBLEM_H
00006
00007 #define DIM 2
00008
00009 #include "LBMProblem.h"
00010 #include "LBMD2Q9Thermal.h"
00011
00012 typedef double DataType;
00013 typedef LBMD2Q9Thermal<DataType> LBMType;
00014
00015 #define OWN_LBMSCHEME
00016 #define OWN_AMRSOLVER
00017 #include "LBMStdProblem.h"
00018 #include "AMRGFMSolver.h"
00019 #include "LBMGFMBoundary.h"
00020 #include "AMRGFMInterpolation.h"
00021 #include "Interfaces/SchemeGFMFileOutput.h"
00022
00023 DataType rd, xs[2], TempS, vr;
00024 DataType pi=4.0*std::atan(1.0);
00025
00026 class LBMSpecific : public LBMType {
00027 typedef LBMType base;
00028 public:
00029 LBMSpecific() : base(), Tmp(1.), deltaTp(1.), gp(1.), betap(1.), Diffp(1.),GasRho(1.), Gasnu(1.), Gasp0(1.) {}
00030
00031 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00032 base::register_at(Ctrl,prefix);
00033 RegisterAt(base::LocCtrl,"Tm",Tmp);
00034 RegisterAt(base::LocCtrl,"deltaT",deltaTp);
00035 RegisterAt(base::LocCtrl,"g",gp);
00036 RegisterAt(base::LocCtrl,"beta",betap);
00037 RegisterAt(base::LocCtrl,"Diff",Diffp);
00038 RegisterAt(base::LocCtrl,"Gas_rho",GasRho);
00039 RegisterAt(base::LocCtrl,"Gas_nu",Gasnu);
00040 RegisterAt(base::LocCtrl,"Gas_p0",Gasp0);
00041 }
00042
00043 virtual void SetupData(GridHierarchy* gh, const int& ghosts) {
00044 base::SetupData(gh,ghosts);
00045
00046 BBox wb = base::GH().wholebbox();
00047 DCoords lc = base::GH().worldCoords(wb.lower(),wb.stepsize());
00048 DCoords uc = base::GH().worldCoords(wb.upper()+wb.stepsize(),wb.stepsize());
00049
00050
00051
00052 DataType Hp=uc(1)-lc(1);
00053
00054 DataType Speed = base::SpeedUp();
00055
00056 DataType viscp = Gasnu;
00057 DataType Pr = viscp/Diffp;
00058
00059 DataType tScale = base::TimeScale();
00060 DataType lScale = base::LengthScale();
00061 DataType THp = Tmp+deltaTp/DataType(2.);
00062 DataType TCp = Tmp-deltaTp/DataType(2.);
00063
00064 DataType Ra = gp*betap*deltaTp*std::pow(Hp,3)/viscp/Diffp;
00065
00066
00067 DataType GasCsp = base::LatticeSpeedOfSound();
00068
00069 base::SetGas(Gasp0,GasRho,viscp,GasCsp);
00070 base::SetThermalGas(TCp,THp,Diffp,gp,betap);
00071
00072 std::cout << "D2Q19Thermal: Gas_rho=" << base::GasDensity() << " Gas_Cs=" << base::GasSpeedofSound()
00073 << " Gas_nu=" << base::GasViscosity() << " Gas_nut=" << base::GasViscosityT() << std::endl;
00074 std::cout << "\nControl Parameters : " << std::endl;
00075 std::cout << "Ra = " << Ra << " Pr = " << Pr << std::endl;
00076 std::cout << "\nTemperature : " << std::endl;
00077 std::cout << "TH = " << THp << " TC = " << TCp << std::endl;
00078 std::cout << "\nFluid values : " << std::endl;
00079 std::cout << "nu = " << viscp << " D = " << Diffp << " g = " << gp << " betap = " << betap << std::endl;
00080 std::cout << "\nGeometry : " << std::endl;
00081 std::cout << "Channel Height = " << Hp << std::endl;
00082 std::cout << "\nLattice Units : " << std::endl;
00083 std::cout << "(Amroc) Omega = " << base::Omega(base::TimeScale()) << " (Amroc) OmegaT = " << base::OmegaT(base::TimeScale()) << std::endl;
00084 std::cout << "(Amroc) Gas_nu = " << base::LatticeViscosity(base::Omega(base::TimeScale())) << " (Amroc) Gas_nuT = " << base::LatticeViscosityT(base::OmegaT(base::TimeScale())) << std::endl;
00085
00086
00087
00088 std::cout << "\nScaling Factors: " << std::endl;
00089 std::cout << "LengthScale = " << base::LengthScale() << std::endl;
00090 std::cout << "TimeScale = " << base::TimeScale() << std::endl;
00091 std::cout << "VelocityScale = " << base::VelocityScale() << std::endl;
00092 std::cout << "SpeedUp = " << base::SpeedUp() << std::endl;
00093
00094 }
00095
00096 protected:
00097 DataType Tmp, deltaTp, gp, betap, Diffp;
00098 DataType GasRho, Gasnu, Gasp0, Ra;
00099 };
00100
00101 class GFMBoundarySpecific : public LBMGFMBoundary<LBMType,DIM> {
00102 typedef LBMGFMBoundary<LBMType,DIM> base;
00103 public:
00104 GFMBoundarySpecific(LBMType &lbm) : base(lbm) {}
00105
00106 virtual void SetBndry(vec_grid_data_type& gdu, const int& nc, const int* idx,
00107 const MicroType* u, const point_type* xc, const DataType* distance,
00108 const point_type* normal, DataType* aux, const int naux, const double t)
00109 { LBM().GFMBCStandard(gdu,base::Type,nc,idx,u,xc,distance,normal,aux,naux,LBMType::Lattice); }
00110
00111 virtual void SetBndryAux(vec_grid_data_type& gdu, grid_data_type& gdphi,
00112 const MicroType* u, DataType* aux, const int& Level,
00113 double t, const int& nc, const int* idx,
00114 const point_type* xc, const DataType* distance,
00115 const point_type* normal) {
00116 DataType xsc[2];
00117 xsc[0] = xs[0]; xsc[1] = xs[1];
00118 for (register int i=0; i<nc; i++) {
00119 DataType alpha = std::atan((xc[i](1)-xsc[1])/(xc[i](0)-xsc[0]));
00120 if (xc[i](0)<xsc[0]) alpha += pi;
00121 aux[i*base::NAux()] = 0.;
00122 aux[i*base::NAux()+1] = 0.;
00123 aux[i*base::NAux()+2] = TempS;
00124 if (vr!=0.) {
00125 aux[i*base::NAux()] = -vr*std::sin(alpha);
00126 aux[i*base::NAux()+1] = vr*std::cos(alpha);
00127 }
00128 }
00129 }
00130 };
00131
00132
00133 class GFMLevelSetSpecific : public GFMLevelSet<DataType,DIM> {
00134 typedef GFMLevelSet<DataType,DIM> base;
00135
00136 public:
00137 typedef base::grid_fct_type grid_fct_type;
00138 typedef base::grid_data_type grid_data_type;
00139
00140 GFMLevelSetSpecific() {}
00141 virtual ~GFMLevelSetSpecific() {}
00142
00143 virtual void SetGrid(grid_data_type& gdphi, const int& Level, const double& t) {
00144 int Nx = gdphi.extents(0), Ny = gdphi.extents(1);
00145 DCoords lb = base::GH().worldCoords(gdphi.lower(), gdphi.stepsize());
00146 DCoords dx = base::GH().worldStep(gdphi.stepsize());
00147 DataType *phi = (DataType *)gdphi.databuffer();
00148
00149 DataType xsc[2];
00150 xsc[0] = xs[0]; xsc[1] = xs[1];
00151 for (register int j=0; j<Ny; j++) {
00152 DataType yd = (DataType(j)+static_cast<DataType>(0.5))*dx(1) + lb(1) - xsc[1];
00153 for (register int i=0; i<Nx; i++) {
00154 DataType xd = (DataType(i)+static_cast<DataType>(0.5))*dx(0) + lb(0) - xsc[0];
00155 phi[i+Nx*j] = std::sqrt(xd*xd+yd*yd) - rd;
00156 }
00157 }
00158 }
00159 };
00160
00161
00162 class GhostFluidMethodSpecific : public GhostFluidMethod<MicroType,DIM> {
00163 typedef GhostFluidMethod<MicroType,DIM> base;
00164 public:
00165 GhostFluidMethodSpecific(boundary_type* boundary, levelset_type* levelset) :
00166 base(boundary,levelset) {
00167 xs[0]=0.; xs[1]=0.; rd=1.; vr=0.; TempS=0.;
00168 }
00169
00170 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00171 base::register_at(Ctrl,prefix);
00172 RegisterAt(base::LocCtrl,"Center(1)",xs[0]);
00173 RegisterAt(base::LocCtrl,"Center(2)",xs[1]);
00174 RegisterAt(base::LocCtrl,"RotationSpeed",vr);
00175 RegisterAt(base::LocCtrl,"RadiusSphere",rd);
00176 RegisterAt(base::LocCtrl,"TemperatureSphere",TempS);
00177 }
00178 };
00179
00180
00181 class SolverSpecific :
00182 public AMRGFMSolver<MicroType,FixupType,FlagType,DIM> {
00183 typedef AMRGFMSolver<MicroType,FixupType,FlagType,DIM> base;
00184 typedef AMRGFMInterpolation<MicroType,FixupType,FlagType,DIM> interpolation_type;
00185 typedef MicroType::InternalDataType DataType;
00186 typedef SchemeGFMFileOutput<LBMType,FixupType,FlagType,DIM> output_type;
00187 typedef interpolation_type::point_type point_type;
00188 public:
00189 SolverSpecific(IntegratorSpecific& integ,
00190 base::initial_condition_type& init,
00191 base::boundary_conditions_type& bc) : base(integ, init, bc) {
00192 SetLevelTransfer(new LBMF77LevelTransfer<LBMType,DIM>(integ.Scheme(), f_prolong, f_restrict));
00193 SetFileOutput(new SchemeGFMFileOutput<LBMType,FixupType,FlagType,DIM>(*this,integ.Scheme()));
00194 SetFixup(new FixupSpecific(integ.Scheme()));
00195 SetFlagging(new FlaggingSpecific(*this,integ.Scheme()));
00196 AddGFM(new GhostFluidMethodSpecific(new GFMBoundarySpecific(integ.LBM()),
00197 new GFMLevelSetSpecific()));
00198 _Interpolation = new interpolation_type(*this);
00199 IntegrateEvery = 1500; Np=50; NPast=500; Uin=0.01;
00200
00201 }
00202
00203 ~SolverSpecific() {
00204 DeleteGFM(_GFM[0]);
00205 delete _LevelTransfer;
00206 delete _Flagging;
00207 delete _Fixup;
00208 delete _FileOutput;
00209 delete _Interpolation;
00210
00211 }
00212
00213 virtual void SetupData() {
00214 base::SetupData();
00215 base::NAMRTimeSteps = 1;
00216 base::AdaptBndTimeInterpolate = 0;
00217 base::Step[0].LastTime = 1.e37;
00218 base::Step[0].VariableTimeStepping = -1;
00219 base::Step[0].dtv[0] = ((LBMIntegrator<LBMType,DIM> &)Integrator_()).LBM().TimeScale();
00220 _Interpolation->SetupData(base::PGH(), base::NGhosts());
00221 }
00222
00223
00224 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00225 base::register_at(Ctrl,prefix);
00226 RegisterAt(base::LocCtrl,"OutputDragEvery",IntegrateEvery);
00227 RegisterAt(base::LocCtrl,"Np",Np);
00228 RegisterAt(base::LocCtrl,"NPast",NPast);
00229 RegisterAt(base::LocCtrl,"StreamVelocity",Uin);
00230 }
00231 virtual void register_at(ControlDevice& Ctrl) {
00232 base::register_at(Ctrl);
00233 }
00234
00235
00236
00237
00238 virtual double Tick(int VariableTimeStepping, const double dtv[], const double cflv[],
00239 int& Rejections) {
00240
00241 double cfl = base::Tick(VariableTimeStepping, dtv, cflv, Rejections);
00242
00243
00244
00245
00246 if (StepsTaken(base::GH(),0)%IntegrateEvery == 0) {
00247
00248 int me = MY_PROC;
00249 double pi = 4.0*std::atan(1.0);
00250
00251 for (register int l=0; l<=FineLevel(base::GH()); l++) {
00252 int Time = CurrentTime(base::GH(),l);
00253 int press_cnt = base::Dim()+4;
00254 ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time, l,
00255 press_cnt, base::t[l]);
00256 }
00257
00258 int Npoints = 4*Np;
00259 DataType* p = new DataType[Npoints];
00260 point_type* xc = new point_type[Npoints];
00261 point_type* normals = new point_type[Npoints];
00262 DataType* dist = (DataType*) 0;
00263 double dalpha = pi/Np;
00264 point_type sum=0.0;
00265
00266 register int n;
00267
00268 for (n=0; n<Npoints; n++) {
00269 double alpha = dalpha*((n-Npoints/2.)+0.5);
00270 xc[n](0) = xs[0]+rd*std::cos(alpha);
00271 xc[n](1) = xs[1]+rd*std::sin(alpha);
00272 }
00273
00274
00275 _Interpolation->PointsValues(base::Work(),base::t[0],Npoints,xc,p);
00276 _Interpolation->ArrayCombine(VizServer,Npoints,p);
00277
00278
00279 _Interpolation->PointsGeom(base::GFM(0).Phi(),base::t[0],Npoints,xc,normals,dist);
00280 _Interpolation->ArrayCombine(VizServer,base::Dim()*Npoints,normals[0].data());
00281
00282
00283
00284
00285 if (me == VizServer) {
00286 double ds = 2.*pi*rd/Npoints;
00287 for (n=0; n<Npoints; n++)
00288 sum += p[n]*normals[n]*ds;
00289 }
00290
00291 delete [] p;
00292 delete [] xc;
00293 delete [] normals;
00294
00295
00296 if (me == VizServer) {
00297 std::ofstream outfile;
00298 std::ostream* out;
00299 std::string fname = "drag.txt";
00300 outfile.open(fname.c_str(), std::ios::out | std::ios::app);
00301 out = new std::ostream(outfile.rdbuf());
00302 *out << base::t[0] << " " << sum(0) << " " << sum(1) << std::endl;
00303 outfile.close();
00304 delete out;
00305 }
00306
00307
00308
00309
00310
00311
00312 for (register int l=0; l<=FineLevel(base::GH()); l++) {
00313 int Time = CurrentTime(base::GH(),l);
00314 int Temp_cnt = base::Dim()+3;
00315 ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time, l,
00316 Temp_cnt, base::t[l]);
00317 }
00318
00319 int NTpoints;
00320 int finelevel = FineLevel(base::GH());
00321
00322 NTpoints = Npoints;
00323 DataType* T = new DataType[NTpoints];
00324 DataType* T2 = new DataType[NTpoints];
00325 DataType* oT = new DataType[NTpoints];
00326 point_type* xc2 = new point_type[NTpoints];
00327 point_type* xc3 = new point_type[NTpoints];
00328 point_type* normals2 = new point_type[NTpoints];
00329 point_type* normals3 = new point_type[NTpoints];
00330 DataType* dist2 = (DataType*) 0;
00331 DataType* dist3 = (DataType*) 0;
00332 double dalphaT = pi/(NTpoints-1.);
00333 double dX = DeltaX(base::GH(),0,finelevel);
00334 double dY = DeltaX(base::GH(),1,finelevel);
00335 double distT = std::sqrt(2.*dX*dY);
00336 double dn = 0.;
00337
00338
00339
00340
00341 for (n=0; n<NTpoints; n++) {
00342 double alphaT = -(pi/2.)+n*dalphaT;
00343 xc2[n](0) = xs[0]+(rd-dn)*std::cos(alphaT);
00344 xc2[n](1) = xs[1]+(rd-dn)*std::sin(alphaT);
00345 }
00346
00347
00348
00349 _Interpolation->PointsValues(base::Work(),base::t[0],NTpoints,xc2,T);
00350 _Interpolation->ArrayCombine(VizServer,NTpoints,T);
00351
00352
00353
00354 _Interpolation->PointsGeom(base::GFM(0).Phi(),base::t[0],NTpoints,xc2,normals2,dist2);
00355 _Interpolation->ArrayCombine(VizServer,base::Dim()*NTpoints,normals2[0].data());
00356
00357
00358
00359
00360 for (n=0; n<NTpoints; n++) {
00361 xc3[n](0) = xc2[n](0)-(distT+dn)*normals2[n](0);
00362 xc3[n](1) = xc2[n](1)-(distT+dn)*normals2[n](1);
00363 }
00364
00365 _Interpolation->PointsValues(base::Work(),base::t[0],NTpoints,xc3,T2);
00366 _Interpolation->ArrayCombine(VizServer,NTpoints,T2);
00367 _Interpolation->PointsGeom(base::GFM(0).Phi(),base::t[0],NTpoints,xc3,normals3,dist3);
00368 _Interpolation->ArrayCombine(VizServer,base::Dim()*NTpoints,normals3[0].data());
00369
00370
00371 if (me == VizServer) {
00372 std::ofstream outfile2;
00373 std::ostream* outTemp;
00374 std::string fname = "nusselt.txt";
00375 outfile2.open(fname.c_str(), std::ios::out | std::ios::app);
00376 outTemp = new std::ostream(outfile2.rdbuf());
00377 *outTemp << base::t[0] << " ";
00378 for (n=0; n<NTpoints; n++){
00379 oT[n] = T2[n]-T[n];
00380 *outTemp << xc2[n](0) << " " << xc2[n](1) << " " << normals2[n](0) << " " << normals2[n](1) << " " << xc3[n](0) << " " << xc3[n](1) << " " << T[n] << " " << T2[n] << " " << oT[n] << " ";
00381 }
00382 *outTemp << std::endl;
00383 outfile2.close();
00384 delete outTemp;
00385 }
00386
00387 delete [] T;
00388 delete [] T2;
00389 delete [] oT;
00390 delete [] xc2;
00391 delete [] normals2;
00392 delete [] xc3;
00393 delete [] normals3;
00394
00395
00396 register int k;
00397 int NPastp = NPast;
00398 DataType* TPast = new DataType[NPastp];
00399 point_type* xPast = new point_type[NPastp];
00400
00401
00402 for (k=0; k<NPastp; k++) {
00403 xPast[k](0) = rd+k*(5*rd/(NPastp-1));
00404 xPast[k](1) = 0.;
00405 }
00406 _Interpolation->PointsValues(base::Work(),base::t[0],NPastp,xPast,TPast);
00407 _Interpolation->ArrayCombine(VizServer,NPastp,TPast);
00408
00409
00410
00411
00412
00413
00414 for (register int l=0; l<=FineLevel(base::GH()); l++) {
00415 int Time = CurrentTime(base::GH(),l);
00416 int uvelPast_cnt = base::Dim();
00417 ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time, l,
00418 uvelPast_cnt, base::t[l]);
00419 }
00420 DataType* uPast = new DataType[NPastp];
00421
00422
00423 _Interpolation->PointsValues(base::Work(),base::t[0],NPastp,xPast,uPast);
00424 _Interpolation->ArrayCombine(VizServer,NPastp,uPast);
00425
00426
00427
00428
00429
00430
00431 for (register int l=0; l<=FineLevel(base::GH()); l++) {
00432 int Time = CurrentTime(base::GH(),l);
00433 int vvelPast_cnt = base::Dim()+1;
00434 ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time, l,
00435 vvelPast_cnt, base::t[l]);
00436 }
00437 DataType* vPast = new DataType[NPastp];
00438
00439
00440 _Interpolation->PointsValues(base::Work(),base::t[0],NPastp,xPast,vPast);
00441 _Interpolation->ArrayCombine(VizServer,NPastp,vPast);
00442
00443
00444 double bottom = (1.*rd/Uin);
00445 double top = (8.*rd/Uin);
00446 if ((base::t[0] >= bottom) && (base::t[0] <= top)){
00447 if (me == VizServer) {
00448 std::ofstream outfile3;
00449 std::ostream* outPast;
00450 std::string fname = "interpolPast.txt";
00451 outfile3.open(fname.c_str(), std::ios::out | std::ios::app);
00452 outPast = new std::ostream(outfile3.rdbuf());
00453 *outPast << base::t[0] << " ";
00454 for (k=0; k<NPastp; k++){
00455 *outPast << xPast[k](0) << " " << xPast[k](1) << " " << uPast[k] << " " << vPast[k] << " " << TPast[k] << " " ;
00456 }
00457 *outPast << std::endl;
00458 outfile3.close();
00459 delete outPast;
00460 }
00461 }
00462
00463 delete [] xPast;
00464 delete [] uPast;
00465 delete [] vPast;
00466 delete [] TPast;
00467
00468 point_type* xPastFFT = new point_type[1];
00469 DataType* vPastFFT = new DataType[1];
00470 xPastFFT[0](0) = 9*rd;
00471 xPastFFT[0](1) = 0.;
00472 _Interpolation->PointsValues(base::Work(),base::t[0],1,xPastFFT,vPastFFT);
00473 _Interpolation->ArrayCombine(VizServer,1,vPastFFT);
00474
00475 if (me == VizServer) {
00476 std::ofstream outfile4;
00477 std::ostream* outFFT;
00478 std::string fname = "FFTPast.txt";
00479 outfile4.open(fname.c_str(), std::ios::out | std::ios::app);
00480 outFFT = new std::ostream(outfile4.rdbuf());
00481 *outFFT << base::t[0] << " " << vPastFFT[0] << std::endl;
00482 outfile4.close();
00483 delete outFFT;
00484 }
00485
00486 delete [] xPastFFT;
00487 delete [] vPastFFT;
00488
00489 }
00490 return cfl;
00491
00492 }
00493
00494 protected:
00495 int IntegrateEvery, Np, NPast;
00496 double Uin;
00497 interpolation_type* _Interpolation;
00498