00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #ifndef AMROC_PROBLEM_H
00012 #define AMROC_PROBLEM_H
00013
00014 #define DIM 2
00015
00016 #include "LBMProblem.h"
00017 #include "LBMD2Q9Thermal_BC.h"
00018
00019 typedef double DataType;
00020 typedef LBMD2Q9Thermal<DataType> LBMType;
00021
00022 #define OWN_LBMSCHEME
00023 #define OWN_AMRSOLVER
00024 #include "LBMStdProblem.h"
00025 #include "AMRGFMSolver.h"
00026
00027 #include "AMRGFMInterpolation.h"
00028 #include "Interfaces/SchemeGFMFileOutput.h"
00029
00030 DataType Kappa;
00031
00032 class LBMSpecific : public LBMType {
00033 typedef LBMType base;
00034 public:
00035 LBMSpecific() : base(), Tmp(1.), deltaTp(1.), gp(1.), betap(1.), Diffp(1.),GasRho(1.), Gasnu(1.), Gasp0(1.), Ra(1.) {}
00036
00037 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00038 base::register_at(Ctrl,prefix);
00039 RegisterAt(base::LocCtrl,"Tm",Tmp);
00040 RegisterAt(base::LocCtrl,"deltaT",deltaTp);
00041 RegisterAt(base::LocCtrl,"g",gp);
00042 RegisterAt(base::LocCtrl,"beta",betap);
00043 RegisterAt(base::LocCtrl,"Diff",Diffp);
00044 RegisterAt(base::LocCtrl,"Gas_rho",GasRho);
00045 RegisterAt(base::LocCtrl,"Gas_nu",Gasnu);
00046 RegisterAt(base::LocCtrl,"Gas_p0",Gasp0);
00047 RegisterAt(base::LocCtrl,"Ra",Ra);
00048 }
00049
00050 virtual void SetupData(GridHierarchy* gh, const int& ghosts) {
00051 base::SetupData(gh,ghosts);
00052
00053 BBox wb = base::GH().wholebbox();
00054 DCoords lc = base::GH().worldCoords(wb.lower(),wb.stepsize());
00055 DCoords uc = base::GH().worldCoords(wb.upper()+wb.stepsize(),wb.stepsize());
00056
00057
00058
00059 DataType Hp=uc(1)-lc(1);
00060
00061 DataType Speed = base::SpeedUp();
00062
00063 DataType viscp = Gasnu;
00064 DataType Pr = viscp/Diffp;
00065
00066 DataType tScale = base::TimeScale();
00067 DataType lScale = base::LengthScale();
00068 DataType THp = Tmp+deltaTp/DataType(2.);
00069 DataType TCp = Tmp-deltaTp/DataType(2.);
00070
00071
00072
00073
00074 DataType GasCsp = base::LatticeSpeedOfSound();
00075
00076 base::SetGas(Gasp0,GasRho,viscp,GasCsp);
00077
00078 base::SetThermalGas(TCp,THp,Diffp,gp,betap);
00079
00080 std::cout << "D2Q19Thermal: Gas_rho=" << base::GasDensity() << " Gas_Cs=" << base::GasSpeedofSound()
00081 << " Gas_nu=" << base::GasViscosity() << " Gas_nut=" << base::GasViscosityT() << std::endl;
00082 std::cout << "\nControl Parameters : " << std::endl;
00083 std::cout << "Ra = " << Ra << " Pr = " << Pr << std::endl;
00084 std::cout << "\nTemperature : " << std::endl;
00085 std::cout << "TH = " << THp << " TC = " << TCp << std::endl;
00086 std::cout << "\nFluid values : " << std::endl;
00087 std::cout << "visc = " << viscp << " Diff = " << Diffp << " g = " << gp << " betap = " << betap << std::endl;
00088 std::cout << "\nGeometry : " << std::endl;
00089 std::cout << "Channel Width = " << Hp << std::endl;
00090 std::cout << "\nLattice Units : " << std::endl;
00091 std::cout << "(Amroc) Omega = " << base::Omega(base::TimeScale()) << " (Amroc) OmegaT = " << base::OmegaT(base::TimeScale()) << std::endl;
00092 std::cout << "(Amroc) Gas_nu = " << base::LatticeViscosity(base::Omega(tScale)) << " (Amroc) Gas_nuT = " << base::LatticeViscosityT(base::OmegaT(tScale)) << std::endl;
00093
00094 std::cout << "Ra Lat = " << base::gbeta*(Hp/lScale)*(Hp/lScale)*(Hp/lScale)/base::LatticeViscosity(base::Omega(tScale))/base::LatticeViscosityT(base::OmegaT(tScale)) << " Ra Phy = " << gp*betap*Hp*Hp*Hp/viscp/Diffp << std::endl;
00095
00096 std::cout << "\nScaling Factors: " << std::endl;
00097 std::cout << "LengthScale = " << base::LengthScale() << std::endl;
00098 std::cout << "TimeScale = " << base::TimeScale() << std::endl;
00099 std::cout << "VelocityScale = " << base::VelocityScale() << std::endl;
00100 std::cout << "SpeedUp = " << base::SpeedUp() << std::endl;
00101
00102 }
00103
00104
00105 protected:
00106 DataType Tmp, deltaTp, gp, betap, Diffp;
00107 DataType GasRho, Gasnu, Gasp0, Ra;
00108 };
00109
00110
00111 class SolverSpecific :
00112 public AMRGFMSolver<MicroType,FixupType,FlagType,DIM> {
00113 typedef AMRGFMSolver<MicroType,FixupType,FlagType,DIM> base;
00114 typedef AMRGFMInterpolation<MicroType,FixupType,FlagType,DIM> interpolation_type;
00115 typedef MicroType::InternalDataType DataType;
00116 typedef SchemeGFMFileOutput<LBMType,FixupType,FlagType,DIM> output_type;
00117 typedef interpolation_type::point_type point_type;
00118 public:
00119 SolverSpecific(IntegratorSpecific& integ,
00120 base::initial_condition_type& init,
00121 base::boundary_conditions_type& bc) : base(integ, init, bc) {
00122 SetLevelTransfer(new LBMF77LevelTransfer<LBMType,DIM>(integ.Scheme(), f_prolong, f_restrict));
00123 SetFileOutput(new SchemeGFMFileOutput<LBMType,FixupType,FlagType,DIM>(*this,integ.Scheme()));
00124 SetFixup(new FixupSpecific(integ.Scheme()));
00125 SetFlagging(new FlaggingSpecific(*this,integ.Scheme()));
00126 _Interpolation = new interpolation_type(*this);
00127 IntegrateEvery = 100; H=100; Tref=0.5, THot=1.0, TCold=0.0;
00128 }
00129
00130 ~SolverSpecific() {
00131 DeleteGFM(_GFM[0]);
00132 delete _LevelTransfer;
00133 delete _Flagging;
00134 delete _Fixup;
00135 delete _FileOutput;
00136 delete _Interpolation;
00137
00138 }
00139
00140 virtual void SetupData() {
00141 base::SetupData();
00142 base::NAMRTimeSteps = 1;
00143 base::AdaptBndTimeInterpolate = 0;
00144 base::Step[0].LastTime = 1.e37;
00145 base::Step[0].VariableTimeStepping = -1;
00146 base::Step[0].dtv[0] = ((LBMIntegrator<LBMType,DIM> &)Integrator_()).LBM().TimeScale();
00147 _Interpolation->SetupData(base::PGH(), base::NGhosts());
00148
00149 Kappa = ((LBMType &)((LBMIntegrator<LBMType,DIM> &)Integrator_()).LBM()).GasViscosityT();
00150 }
00151
00152
00153 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00154 base::register_at(Ctrl,prefix);
00155 RegisterAt(base::LocCtrl,"OutputEvery",IntegrateEvery);
00156 RegisterAt(base::LocCtrl,"H",H);
00157 RegisterAt(base::LocCtrl,"Tref",Tref);
00158 RegisterAt(base::LocCtrl,"TH",THot);
00159 RegisterAt(base::LocCtrl,"TC",TCold);
00160 }
00161 virtual void register_at(ControlDevice& Ctrl) {
00162 base::register_at(Ctrl);
00163 }
00164
00165 virtual double Tick(int VariableTimeStepping, const double dtv[], const double cflv[],
00166 int& Rejections) {
00167
00168 double cfl = base::Tick(VariableTimeStepping, dtv, cflv, Rejections);
00169
00170 if (StepsTaken(base::GH(),0)%IntegrateEvery == 0) {
00171
00172 int me = MY_PROC;
00173
00174
00175
00176
00177
00178
00179 for (register int l=0; l<=FineLevel(base::GH()); l++) {
00180 int Time = CurrentTime(base::GH(),l);
00181 int Temp_cnt = base::Dim()+3;
00182 ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time, l,
00183 Temp_cnt, base::t[l]);
00184 }
00185
00186 double deltaT = THot-TCold;
00187 int NTpoints;
00188 int finelevel = FineLevel(base::GH());
00189
00190 double dX = DeltaX(base::GH(),0,finelevel);
00191 double dY = DeltaX(base::GH(),1,finelevel);
00192 NTpoints = H/dX;
00193
00194 DataType* T = new DataType[NTpoints];
00195 DataType* T2 = new DataType[NTpoints];
00196 DataType* oT = new DataType[NTpoints];
00197 point_type* xc = new point_type[NTpoints];
00198 point_type* xc2 = new point_type[NTpoints];
00199
00200 std::cout << "NTpoints = " << NTpoints << ", fineLevel = " << finelevel << ", dX/dY = " << dX << " " << dY << std::endl;
00201
00202
00203 register int n;
00204
00205 for (n=0; n<NTpoints; n++) {
00206 xc[n](0) = dX/2.;
00207 xc[n](1) = (dY/2.)+n*dY;
00208 xc2[n](0) = xc[n](0)+dX;
00209 xc2[n](1) = xc[n](1);
00210 }
00211
00212
00213
00214 _Interpolation->PointsValues(base::Work(),base::t[0],NTpoints,xc,T);
00215 _Interpolation->ArrayCombine(VizServer,NTpoints,T);
00216
00217 _Interpolation->PointsValues(base::Work(),base::t[0],NTpoints,xc2,T2);
00218 _Interpolation->ArrayCombine(VizServer,NTpoints,T2);
00219
00220
00221 if (me == VizServer) {
00222 std::ofstream outfile;
00223 std::ostream* outTemp;
00224 std::string fname = "TemperatureFull.txt";
00225 outfile.open(fname.c_str(), std::ios::out | std::ios::app);
00226 outTemp = new std::ostream(outfile.rdbuf());
00227 outTemp->precision(12);
00228 *outTemp << base::t[0] << " ";
00229 for (n=0; n<NTpoints; n++){
00230 oT[n] = (T2[n]-T[n])/(xc2[n](0)-xc[n](0));
00231 *outTemp << xc[n](0) << " " << xc[n](1) << " " << T[n] << " " << xc2[n](0) << " " << xc2[n](1) << " " << T2[n] << " " << oT[n] << " ";
00232 }
00233 *outTemp << std::endl;
00234 outfile.close();
00235 delete outTemp;
00236
00237 DataType nusselt, nusseltL;
00238 nusselt = DataType(0.);
00239 nusseltL = DataType(0.);
00240
00241 std::ofstream outfile3;
00242 std::ostream* outTemp3;
00243 std::string fname3 = "Nusselt.txt";
00244 outfile3.open(fname3.c_str(), std::ios::out | std::ios::app);
00245 outTemp3 = new std::ostream(outfile3.rdbuf());
00246 outTemp3->precision(12);
00247 *outTemp3 << base::t[0] << " ";
00248 for (n=0; n<NTpoints; n++){
00249 nusseltL = oT[n];
00250 nusselt += nusseltL;
00251 }
00252 nusselt = -nusselt/(NTpoints-1)*(xc[NTpoints-1](1)-xc[0](1));
00253 *outTemp3 << nusselt << std::endl;
00254 outfile3.close();
00255 delete outTemp3;
00256
00257 }
00258
00259 delete [] T;
00260 delete [] T2;
00261 delete [] oT;
00262 delete [] xc2;
00263
00264
00265 DataType* TW = new DataType[NTpoints];
00266 DataType* TA = new DataType[NTpoints];
00267 DataType* TB = new DataType[NTpoints];
00268 DataType* NussExtr = new DataType[NTpoints];
00269 point_type* xW = new point_type[NTpoints];
00270 point_type* xA = new point_type[NTpoints];
00271 point_type* xB = new point_type[NTpoints];
00272
00273
00274 for (n=0; n<NTpoints; n++) {
00275 xW[n](0) = 0.;
00276 xW[n](1) = (dY/2.)+n*dY;
00277 xA[n](0) = xc[n](0);
00278 xA[n](1) = xW[n](1);
00279 xB[n](0) = xc[n](0)+dX;
00280 xB[n](1) = xW[n](1);
00281 }
00282
00283 _Interpolation->PointsValues(base::Work(),base::t[0],NTpoints,xA,TA);
00284 _Interpolation->ArrayCombine(VizServer,NTpoints,TA);
00285 _Interpolation->PointsValues(base::Work(),base::t[0],NTpoints,xB,TB);
00286 _Interpolation->ArrayCombine(VizServer,NTpoints,TB);
00287
00288 if (me == VizServer) {
00289
00290 std::ofstream outfileExt;
00291 std::ostream* outTempExt;
00292 std::string fname = "TempExt.txt";
00293 outfileExt.open(fname.c_str(), std::ios::out | std::ios::app);
00294 outTempExt = new std::ostream(outfileExt.rdbuf());
00295 outTempExt->precision(12);
00296 for (n=0; n<NTpoints; n++){
00297 TW[n] = TB[n]+(xW[n](0)-xB[n](0))/(xA[n](0)-xB[n](0))*(TA[n]-TB[n]);
00298 NussExtr[n] = -(TA[n]-TW[n])/(xA[n](0)-xW[n](0));
00299 *outTempExt << xW[n](0) << " " << xW[n](1) << " " << TW[n] << " " << xA[n](0) << " " << xA[n](1) << " " << TA[n] << " " << xB[n](0) << " " << xB[n](1) << " " << TB[n] << " " << NussExtr[n] << " ";
00300 }
00301 *outTempExt << std::endl;
00302 outfileExt.close();
00303 delete outTempExt;
00304
00305
00306 DataType NussAveExt;
00307
00308 std::ofstream outfileNuss;
00309 std::ostream* outNussExt;
00310 std::string fname4 = "NussExt.txt";
00311 outfileExt.open(fname4.c_str(), std::ios::out | std::ios::app);
00312 outNussExt = new std::ostream(outfileExt.rdbuf());
00313 outNussExt->precision(12);
00314 for (n=0; n<NTpoints; n++){
00315 NussAveExt += NussExtr[n];
00316 }
00317 NussAveExt = NussAveExt/(NTpoints-1)*(xc[NTpoints-1](1)-xc[0](1));
00318 *outNussExt << NussAveExt << " " << std::endl;
00319 outfileNuss.close();
00320 delete outNussExt;
00321 }
00322
00323 delete [] TW;
00324 delete [] TA;
00325 delete [] TB;
00326 delete [] xW;
00327 delete [] xA;
00328 delete [] xB;
00329 delete [] xc;
00330
00331
00332
00333
00334
00335
00336
00337 for(register int l=0; l<=FineLevel(base::GH()); l++) {
00338 int Time = CurrentTime(base::GH(),l);
00339 int Temp_cnt = base::Dim();
00340 ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time, l,
00341 Temp_cnt, base::t[l]);
00342 }
00343
00344 int Nupoints;
00345 Nupoints = H/dX;
00346
00347 DataType* u = new DataType[Nupoints];
00348 DataType* u2 = new DataType[Nupoints];
00349 DataType* um = new DataType[Nupoints];
00350 point_type* yu = new point_type[Nupoints];
00351 point_type* yu2 = new point_type[Nupoints];
00352
00353
00354
00355 for (n=0; n<Nupoints; n++) {
00356 yu[n](1) = dY/2.+n*dY;
00357 yu[n](0) = (H/2)-(dY/2);
00358 yu2[n](1) = yu[n](1);
00359 yu2[n](0) = (H/2)+(dY/2);
00360 }
00361
00362
00363
00364 _Interpolation->PointsValues(base::Work(),base::t[0],Nupoints,yu,u);
00365 _Interpolation->ArrayCombine(VizServer,Nupoints,u);
00366
00367 _Interpolation->PointsValues(base::Work(),base::t[0],Nupoints,yu2,u2);
00368 _Interpolation->ArrayCombine(VizServer,Nupoints,u2);
00369
00370
00371 if (me == VizServer) {
00372
00373 DataType* uMax;
00374 point_type* yMax = new point_type[1];
00375
00376 std::ofstream outfileU;
00377 std::ostream* outU;
00378 std::string fname5 = "Vel_u.txt";
00379 outfileU.open(fname5.c_str(), std::ios::out | std::ios::app);
00380 outU = new std::ostream(outfileU.rdbuf());
00381 outU->precision(12);
00382 for (n=0; n<Nupoints; n++){
00383
00384 um[n] = (u[n]+u2[n])/2.;
00385 *outU << yu[n](0)/H << " " << yu[n](1)/H << " " << u[n]*H/Kappa << " " << yu2[n](0)/H << " " << yu2[n](1)/H << " " << u2[n]*H/Kappa << " " << um[n]*H/Kappa << " ";
00386 }
00387 *outU << std::endl;
00388 outfileU.close();
00389 delete outU;
00390
00391 uMax = std::max_element(um,um+Nupoints);
00392 yMax[0](0) = (yu[uMax-um](0)+yu2[uMax-um](0))/2.;
00393 yMax[0](1) = yu[uMax-um](1);
00394
00395 std::ofstream outfileUmax;
00396 std::ostream* outUmax;
00397 fname5 = "Vel_uMax.txt";
00398 outfileUmax.open(fname5.c_str(), std::ios::out | std::ios::app);
00399 outUmax = new std::ostream(outfileUmax.rdbuf());
00400 outUmax->precision(12);
00401 *outUmax << *uMax*H/Kappa << " " << yMax[0](0)/H << " " << yMax[0](1)/H << std::endl;
00402 outfileUmax.close();
00403 delete outUmax;
00404 }
00405
00406
00407
00408
00409
00410
00411
00412 for(register int l=0; l<=FineLevel(base::GH()); l++) {
00413 int Time = CurrentTime(base::GH(),l);
00414 int Temp_cnt = base::Dim()+1;
00415 ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time, l,
00416 Temp_cnt, base::t[l]);
00417 }
00418
00419 int Nvpoints;
00420 Nvpoints = H/dX;
00421
00422 DataType* v = new DataType[Nvpoints];
00423 DataType* v2 = new DataType[Nvpoints];
00424 DataType* vm = new DataType[Nvpoints];
00425 point_type* xv = new point_type[Nvpoints];
00426 point_type* xv2 = new point_type[Nvpoints];
00427
00428
00429
00430 for (n=0; n<Nvpoints; n++) {
00431 xv[n](0) = dX/2.+n*dX;
00432 xv[n](1) = (H/2)-(dX/2);
00433 xv2[n](0) = xv[n](0);
00434 xv2[n](1) = (H/2)+(dX/2);
00435 }
00436
00437
00438
00439 _Interpolation->PointsValues(base::Work(),base::t[0],Nvpoints,xv,v);
00440 _Interpolation->ArrayCombine(VizServer,Nupoints,v);
00441
00442 _Interpolation->PointsValues(base::Work(),base::t[0],Nvpoints,xv2,v2);
00443 _Interpolation->ArrayCombine(VizServer,Nvpoints,v2);
00444
00445
00446 if (me == VizServer) {
00447
00448 DataType* vMax;
00449 point_type* xMax = new point_type[1];
00450
00451 std::ofstream outfileV;
00452 std::ostream* outV;
00453 std::string fname5 = "Vel_v.txt";
00454 outfileV.open(fname5.c_str(), std::ios::out | std::ios::app);
00455 outV = new std::ostream(outfileV.rdbuf());
00456 outV->precision(12);
00457 for (n=0; n<Nvpoints; n++){
00458 vm[n] = (v[n]+v2[n])/2.;
00459 *outV << xv[n](0)/H << " " << xv[n](1)/H << " " << v[n]*H/Kappa << " " << xv2[n](0)/H << " " << xv2[n](1)/H << " " << v2[n]*H/Kappa << " " << vm[n]*H/Kappa << " ";
00460 }
00461 *outV << std::endl;
00462 outfileV.close();
00463 delete outV;
00464
00465 vMax = std::max_element(vm,vm+Nvpoints);
00466 xMax[0](1) = (xv[vMax-vm](1)+xv2[vMax-vm](1))/2.;
00467 xMax[0](0) = xv[vMax-vm](0);
00468
00469 std::ofstream outfileVmax;
00470 std::ostream* outVmax;
00471 fname5 = "Vel_vMax.txt";
00472 outfileVmax.open(fname5.c_str(), std::ios::out | std::ios::app);
00473 outVmax = new std::ostream(outfileVmax.rdbuf());
00474 outVmax->precision(12);
00475 *outVmax << *vMax*H/Kappa << " " << xMax[0](0)/H << " " << xMax[0](1)/H << std::endl;
00476 outfileVmax.close();
00477 delete outVmax;
00478 }
00479
00480
00481
00482 }
00483
00484
00485 return cfl;
00486
00487 }
00488
00489 protected:
00490 int IntegrateEvery, H;
00491 double Tref, THot, TCold;
00492 interpolation_type* _Interpolation;
00493
00494 };
00495 #endif
00496
00497