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 "Interfaces/SchemeGFMFileOutput.h"
00021
00022 DataType p0 = 1.;
00023 DataType rd, xs[2], rs, T, Tr,TempS;
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.), Ra(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 RegisterAt(base::LocCtrl,"Ra",Ra);
00042 }
00043
00044 virtual void SetupData(GridHierarchy* gh, const int& ghosts) {
00045 base::SetupData(gh,ghosts);
00046
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
00051
00052
00053 DataType Hp=uc(1)-lc(1);
00054
00055 DataType Speed = base::SpeedUp();
00056
00057 DataType viscp = Gasnu;
00058 DataType Pr = viscp/Diffp;
00059
00060 DataType tScale = base::TimeScale();
00061 DataType lScale = base::LengthScale();
00062 DataType THp = Tmp+deltaTp/DataType(2.);
00063 DataType TCp = Tmp-deltaTp/DataType(2.);
00064
00065
00066
00067
00068 DataType GasCsp = base::LatticeSpeedOfSound();
00069
00070 base::SetGas(Gasp0,GasRho,viscp,GasCsp);
00071
00072 base::SetThermalGas(TCp,THp,Diffp,gp,betap);
00073
00074 std::cout << "D2Q19Thermal: Gas_rho=" << base::GasDensity() << " Gas_Cs=" << base::GasSpeedofSound()
00075 << " Gas_nu=" << base::GasViscosity() << " Gas_nut=" << base::GasViscosityT() << std::endl;
00076 std::cout << "\nControl Parameters : " << std::endl;
00077 std::cout << "Ra = " << Ra << " Pr = " << Pr << std::endl;
00078 std::cout << "\nTemperature : " << std::endl;
00079 std::cout << "TH = " << THp << " TC = " << TCp << std::endl;
00080 std::cout << "\nFluid values : " << std::endl;
00081 std::cout << "visc = " << viscp << " Diff = " << Diffp << " g = " << gp << " betap = " << betap << std::endl;
00082 std::cout << "\nGeometry : " << std::endl;
00083 std::cout << "Channel Width = " << Hp << std::endl;
00084 std::cout << "\nLattice Units : " << std::endl;
00085 std::cout << "(Amroc) Omega = " << base::Omega(base::TimeScale()) << " (Amroc) OmegaT = " << base::OmegaT(base::TimeScale()) << std::endl;
00086 std::cout << "(Amroc) Gas_nu = " << base::LatticeViscosity(base::Omega(base::TimeScale())) << " (Amroc) Gas_nuT = " << base::LatticeViscosityT(base::OmegaT(base::TimeScale())) << std::endl;
00087
00088
00089
00090 std::cout << "\nScaling Factors: " << std::endl;
00091 std::cout << "LengthScale = " << base::LengthScale() << std::endl;
00092 std::cout << "TimeScale = " << base::TimeScale() << std::endl;
00093 std::cout << "VelocityScale = " << base::VelocityScale() << std::endl;
00094 std::cout << "SpeedUp = " << base::SpeedUp() << std::endl;
00095
00096 }
00097
00098
00099 protected:
00100 DataType Tmp, deltaTp, gp, betap, Diffp;
00101 DataType GasRho, Gasnu, Gasp0, Ra;
00102 };
00103
00104 class GFMBoundarySpecific : public LBMGFMBoundary<LBMType,DIM> {
00105 typedef LBMGFMBoundary<LBMType,DIM> base;
00106 public:
00107 GFMBoundarySpecific(LBMType &lbm) : base(lbm) {}
00108
00109 virtual void SetBndry(vec_grid_data_type& gdu, const int& nc, const int* idx,
00110 const MicroType* u, const point_type* xc, const DataType* distance,
00111 const point_type* normal, DataType* aux, const int naux, const double t)
00112 { LBM().GFMBCStandard(gdu,base::Type,nc,idx,u,xc,distance,normal,aux,naux,LBMType::Lattice); }
00113
00114 virtual void SetBndryAux(vec_grid_data_type& gdu, grid_data_type& gdphi,
00115 const MicroType* u, DataType* aux, const int& Level,
00116 double t, const int& nc, const int* idx,
00117 const point_type* xc, const DataType* distance,
00118 const point_type* normal) {
00119 DataType xsc[2], vc[2];
00120 xsc[0] = xs[0]; xsc[1] = xs[1];
00121 vc[0] = 0.; vc[1] = 0.;
00122 if (T!=0.) {
00123 xsc[0] += rs*std::cos(2.*pi/T*t);
00124 xsc[1] += rs*std::sin(2.*pi/T*t);
00125 vc[0] = -2.*pi/T*std::sin(2.*pi/T*t);
00126 vc[1] = 2.*pi/T*std::cos(2.*pi/T*t);
00127 }
00128 for (register int i=0; i<nc; i++) {
00129
00130
00131 DataType alpha = std::atan((xc[i](1)-xsc[1])/(xc[i](0)-xsc[0]));
00132 if (xc[i](0)<xsc[0]) alpha += pi;
00133 DataType xp = xsc[0]+rd*std::cos(alpha)-xs[0];
00134 DataType yp = xsc[1]+rd*std::sin(alpha)-xs[1];
00135 DataType rc = std::sqrt(xp*xp+yp*yp);
00136 aux[i*base::NAux()] = rc*vc[0];
00137 aux[i*base::NAux()+1] = rc*vc[1];
00138 aux[i*base::NAux()+2] = TempS;
00139 if (Tr!=0.) {
00140 DataType vr = 2.*pi/Tr*rd;
00141 aux[i*base::NAux()] -= vr*std::sin(alpha);
00142 aux[i*base::NAux()+1] += vr*std::cos(alpha);
00143 aux[i*base::NAux()+2] = TempS;
00144 }
00145 }
00146 }
00147 };
00148
00149
00150 class GFMLevelSetSpecific : public GFMLevelSet<DataType,DIM> {
00151 typedef GFMLevelSet<DataType,DIM> base;
00152
00153 public:
00154 typedef base::grid_fct_type grid_fct_type;
00155 typedef base::grid_data_type grid_data_type;
00156
00157 GFMLevelSetSpecific() {}
00158 virtual ~GFMLevelSetSpecific() {}
00159
00160 virtual void SetGrid(grid_data_type& gdphi, const int& Level, const double& t) {
00161 int Nx = gdphi.extents(0), Ny = gdphi.extents(1);
00162 DCoords lb = base::GH().worldCoords(gdphi.lower(), gdphi.stepsize());
00163 DCoords dx = base::GH().worldStep(gdphi.stepsize());
00164 DataType *phi = (DataType *)gdphi.databuffer();
00165
00166 DataType xsc[2];
00167 xsc[0] = xs[0]; xsc[1] = xs[1];
00168 if (T!=0.) {
00169 xsc[0] += rs*std::cos(2.*pi/T*t);
00170 xsc[1] += rs*std::sin(2.*pi/T*t);
00171 }
00172 for (register int j=0; j<Ny; j++) {
00173 DataType yd = (DataType(j)+static_cast<DataType>(0.5))*dx(1) + lb(1) - xsc[1];
00174 for (register int i=0; i<Nx; i++) {
00175 DataType xd = (DataType(i)+static_cast<DataType>(0.5))*dx(0) + lb(0) - xsc[0];
00176 phi[i+Nx*j] = std::sqrt(xd*xd+yd*yd) - rd;
00177 }
00178 }
00179 }
00180 };
00181
00182
00183 class GhostFluidMethodSpecific : public GhostFluidMethod<MicroType,DIM> {
00184 typedef GhostFluidMethod<MicroType,DIM> base;
00185 public:
00186 GhostFluidMethodSpecific(boundary_type* boundary, levelset_type* levelset) :
00187 base(boundary,levelset) {
00188 xs[0]=0.; xs[1]=0.; rs = 1.; rd=1.; T=100.; Tr=0.; TempS=0.;
00189 }
00190
00191 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00192 base::register_at(Ctrl,prefix);
00193 RegisterAt(base::LocCtrl,"Center(1)",xs[0]);
00194 RegisterAt(base::LocCtrl,"Center(2)",xs[1]);
00195 RegisterAt(base::LocCtrl,"RotationDuration",Tr);
00196 RegisterAt(base::LocCtrl,"CircleDuration",T);
00197 RegisterAt(base::LocCtrl,"RadiusCircle",rs);
00198 RegisterAt(base::LocCtrl,"RadiusSphere",rd);
00199 RegisterAt(base::LocCtrl,"TemperatureSphere",TempS);
00200 }
00201 };
00202
00203
00204 class SolverSpecific :
00205 public AMRGFMSolver<MicroType,FixupType,FlagType,DIM> {
00206 typedef AMRGFMSolver<MicroType,FixupType,FlagType,DIM> base;
00207 public:
00208 SolverSpecific(IntegratorSpecific& integ,
00209 base::initial_condition_type& init,
00210 base::boundary_conditions_type& bc) : base(integ, init, bc) {
00211 SetLevelTransfer(new LBMF77LevelTransfer<LBMType,DIM>(integ.Scheme(), f_prolong, f_restrict));
00212 SetFileOutput(new SchemeGFMFileOutput<LBMType,FixupType,FlagType,DIM>(*this,integ.Scheme()));
00213 SetFixup(new FixupSpecific(integ.Scheme()));
00214 SetFlagging(new FlaggingSpecific(*this,integ.Scheme()));
00215 AddGFM(new GhostFluidMethodSpecific(new GFMBoundarySpecific(integ.LBM()),
00216 new GFMLevelSetSpecific()));
00217 }
00218
00219 ~SolverSpecific() {
00220 DeleteGFM(_GFM[0]);
00221 delete _LevelTransfer;
00222 delete _Flagging;
00223 delete _Fixup;
00224 delete _FileOutput;
00225 }
00226
00227 virtual void SetupData() {
00228 base::SetupData();
00229 base::NAMRTimeSteps = 1;
00230 base::AdaptBndTimeInterpolate = 0;
00231 base::Step[0].LastTime = 1.e37;
00232 base::Step[0].VariableTimeStepping = -1;
00233 base::Step[0].dtv[0] = ((LBMIntegrator<LBMType,DIM> &)Integrator_()).LBM().TimeScale();
00234 }
00235 };