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 #define OWN_BOUNDARYCONDITION
00016 #define OWN_INITIALCONDITION
00017 #include "LBMStdProblem.h"
00018 #include "LBMD3Q19SupplementalBC.h"
00019 #include "AMRGFMSolver.h"
00020 #include "LBMGFMBoundary.h"
00021 #include "AMRGFMInterpolation.h"
00022 #include "Interfaces/SchemeGFMFileOutput.h"
00023
00024 DataType xs[3], rd, u0, r0, nu, p0;
00025
00026 class GFMLevelSetSpecific : public GFMLevelSet<DataType,DIM> {
00027 typedef GFMLevelSet<DataType,DIM> base;
00028
00029 public:
00030 typedef base::grid_fct_type grid_fct_type;
00031 typedef base::grid_data_type grid_data_type;
00032
00033 GFMLevelSetSpecific() {}
00034 virtual ~GFMLevelSetSpecific() {}
00035
00036 virtual void SetGrid(grid_data_type& gdphi, const int& Level, const double& t) {
00037 int Nx = gdphi.extents(0), Ny = gdphi.extents(1), Nz = gdphi.extents(2);
00038 DCoords lb = base::GH().worldCoords(gdphi.lower(), gdphi.stepsize());
00039 DCoords dx = base::GH().worldStep(gdphi.stepsize());
00040 DataType *phi = (DataType *)gdphi.databuffer();
00041
00042 for (register int k=0; k<Nz; k++) {
00043 DataType zd = (DataType(k)+static_cast<DataType>(0.5))*dx(2) + lb(2) - xs[2];
00044 for (register int j=0; j<Ny; j++) {
00045 DataType yd = (DataType(j)+static_cast<DataType>(0.5))*dx(1) + lb(1) - xs[1];
00046 for (register int i=0; i<Nx; i++) {
00047 DataType xd = (DataType(i)+static_cast<DataType>(0.5))*dx(0) + lb(0) - xs[0];
00048 phi[i+Nx*(j+Ny*k)] = std::sqrt(xd*xd+yd*yd+zd*zd) - rd;
00049 }
00050 }
00051 }
00052 }
00053 };
00054
00055
00056 class GhostFluidMethodSpecific : public GhostFluidMethod<MicroType,DIM> {
00057 typedef GhostFluidMethod<MicroType,DIM> base;
00058 public:
00059 GhostFluidMethodSpecific(boundary_type* boundary, levelset_type* levelset) :
00060 base(boundary,levelset) {
00061 xs[0]=0.; xs[1]=0.; xs[2]=0.; rd=0.;
00062 }
00063
00064 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00065 base::register_at(Ctrl,prefix);
00066 RegisterAt(base::LocCtrl,"Center(1)",xs[0]);
00067 RegisterAt(base::LocCtrl,"Center(2)",xs[1]);
00068 RegisterAt(base::LocCtrl,"Center(3)",xs[2]);
00069 RegisterAt(base::LocCtrl,"Radius",rd);
00070 }
00071 };
00072
00073
00074 class BoundaryConditionsSpecific : public BoundaryConditionsSupplemental {
00075 public:
00076 BoundaryConditionsSpecific(LBMType &lbm) : BoundaryConditionsSupplemental(lbm) {}
00077 };
00078
00079
00080 class InitialConditionSpecific : public LBMInitialCondition<LBMType,DIM> {
00081 typedef LBMInitialCondition<LBMType,DIM> base;
00082
00083 public:
00084 InitialConditionSpecific(LBMType &lbm) : base(lbm) {}
00085
00086 virtual void SetupData(GridHierarchy* gh, const int& ghosts) {
00087 base::SetupData(gh,ghosts);
00088 if (base::NAux>=2) {
00089 r0 = aux[0]; u0 = aux[1];
00090 }
00091 }
00092 };
00093
00094
00095 class SolverSpecific :
00096 public AMRGFMSolver<MicroType,FixupType,FlagType,DIM> {
00097 typedef AMRGFMSolver<MicroType,FixupType,FlagType,DIM> base;
00098 typedef AMRGFMInterpolation<MicroType,FixupType,FlagType,DIM> interpolation_type;
00099 typedef MicroType::InternalDataType DataType;
00100 typedef SchemeGFMFileOutput<LBMType,FixupType,FlagType,DIM> output_type;
00101 typedef interpolation_type::point_type point_type;
00102 public:
00103 SolverSpecific(IntegratorSpecific& integ,
00104 InitialConditionSpecific& init,
00105 base::boundary_conditions_type& bc) : base(integ, init, bc) {
00106 SetLevelTransfer(new LBMF77LevelTransfer<LBMType,DIM>(integ.Scheme(), f_prolong, f_restrict));
00107 SetFileOutput(new SchemeGFMFileOutput<LBMType,FixupType,FlagType,DIM>(*this,integ.Scheme()));
00108 SetFixup(new FixupSpecific(integ.Scheme()));
00109 SetFlagging(new FlaggingSpecific(*this,integ.Scheme()));
00110 AddGFM(new GhostFluidMethodSpecific(new LBMGFMBoundary<LBMType,DIM>(integ.LBM()),
00111 new GFMLevelSetSpecific()));
00112 _Interpolation = new interpolation_type(*this);
00113 IntegrateEvery = 1; Np=100;
00114 SurfaceFile = "-";
00115 pi = 4.0*std::atan(1.0);
00116 dev = (DataType *)0; sig = (DataType *)0;
00117 p = (DataType *)0; pa = (DataType *)0;
00118 taut = (DataType *)0; tauta = (DataType *)0;
00119 tauv = (DataType *)0; tauva = (DataType *)0;
00120 xc = (point_type *)0;
00121 pac = 0;
00122 }
00123
00124 ~SolverSpecific() {
00125 DeleteGFM(_GFM[0]);
00126 delete _LevelTransfer;
00127 delete _Flagging;
00128 delete _Fixup;
00129 delete _FileOutput;
00130 delete _Interpolation;
00131 if (dev) delete [] dev;
00132 if (sig) delete [] sig;
00133 if (p) delete [] p;
00134 if (pa) delete [] pa;
00135 if (taut) delete [] taut;
00136 if (tauta) delete [] tauta;
00137 if (tauv) delete [] tauv;
00138 if (tauva) delete [] tauva;
00139 if (xc) delete [] xc;
00140 }
00141
00142 virtual void SetupData() {
00143 base::SetupData();
00144 base::NAMRTimeSteps = 1;
00145 base::AdaptBndTimeInterpolate = 0;
00146 base::Step[0].LastTime = 1.e37;
00147 base::Step[0].VariableTimeStepping = -1;
00148 base::Step[0].dtv[0] = ((LBMIntegrator<LBMType,DIM> &)Integrator_()).LBM().TimeScale();
00149 _Interpolation->SetupData(base::PGH(), base::NGhosts());
00150
00151 int Npoints = 8*Np*Np;
00152 dev = new DataType[6*Npoints];
00153 sig = new DataType[6*Npoints];
00154 p = new DataType[Npoints]; pa = new DataType[Npoints];
00155 taut = new DataType[Npoints]; tauta = new DataType[Npoints];
00156 tauv = new DataType[Npoints]; tauva = new DataType[Npoints];
00157 for (register int i=0; i<Npoints; i++) {
00158 pa[i]=0.; tauta[i]=0.; tauva[i]=0.;
00159 }
00160 xc = new point_type[Npoints];
00161 nu = ((LBMType &)((LBMIntegrator<LBMType,DIM> &)Integrator_()).LBM()).GasViscosity();
00162 p0 = ((LBMType &)((LBMIntegrator<LBMType,DIM> &)Integrator_()).LBM()).BasePressure();
00163
00164 std::cout << "r0=" << r0 << " u0=" << u0 << " rd=" << rd << " p0=" << p0
00165 << " Re=" << u0*2.*rd/nu << std::endl;
00166 }
00167
00168 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00169 base::register_at(Ctrl,prefix);
00170 RegisterAt(base::LocCtrl,"OutputDragEvery",IntegrateEvery);
00171 RegisterAt(base::LocCtrl,"Np",Np);
00172 RegisterAt(base::LocCtrl,"SurfacePlot",SurfaceFile);
00173 }
00174 virtual void register_at(ControlDevice& Ctrl) {
00175 base::register_at(Ctrl);
00176 }
00177
00178 void EvalForce(point_type& st, point_type& sp, point_type& sv) {
00179 int me = MY_PROC;
00180 int Npoints = 8*Np*Np;
00181 point_type* normals = new point_type[Npoints];
00182 DataType* ds = new DataType[Npoints];
00183 DataType* dist = (DataType*) 0;
00184 double dalpha = 0.5*pi/Np;
00185 double dbeta = 0.5*pi/Np;
00186 register int n = 0;
00187
00188
00189 for (register int l=0; l<4*Np; l++) {
00190 double alpha = dalpha*(l+0.5);
00191 for (register int m=0; m<2*Np; m++) {
00192 double beta = dalpha*(m+0.5);
00193 xc[n](0) = xs[0]+rd*std::cos(alpha)*std::sin(beta);
00194 xc[n](1) = xs[1]+rd*std::sin(alpha)*std::sin(beta);
00195 xc[n](2) = xs[2]+rd*std::cos(beta);
00196 ds[n] = pi*rd*rd*dbeta*std::sin(beta)/(2.*Np);
00197 n++;
00198 }
00199 }
00200
00201
00202 for (register int l=0; l<=FineLevel(base::GH()); l++) {
00203 int Time = CurrentTime(base::GH(),l);
00204 ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time, l,
00205 7, base::t[l]);
00206 }
00207
00208 _Interpolation->PointsValues(base::Work(),base::t[0],Npoints,xc,p);
00209 _Interpolation->ArrayCombine(VizServer,Npoints,p);
00210
00211
00212 for (register int cnt=19; cnt<=24; cnt++) {
00213 for (register int l=0; l<=FineLevel(base::GH()); l++) {
00214 int Time = CurrentTime(base::GH(),l);
00215 ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time, l,
00216 cnt, base::t[l]);
00217 }
00218
00219 _Interpolation->PointsValues(base::Work(),base::t[0],Npoints,xc,dev+(cnt-19)*Npoints);
00220 }
00221 _Interpolation->ArrayCombine(VizServer,6*Npoints,dev);
00222
00223
00224 for (register int cnt=28; cnt<=33; cnt++) {
00225 for (register int l=0; l<=FineLevel(base::GH()); l++) {
00226 int Time = CurrentTime(base::GH(),l);
00227 ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time, l,
00228 cnt, base::t[l]);
00229 }
00230
00231 _Interpolation->PointsValues(base::Work(),base::t[0],Npoints,xc,sig+(cnt-28)*Npoints);
00232 }
00233 _Interpolation->ArrayCombine(VizServer,6*Npoints,sig);
00234
00235
00236 _Interpolation->PointsGeom(base::GFM(0).Phi(),base::t[0],Npoints,xc,normals,dist);
00237 _Interpolation->ArrayCombine(VizServer,base::Dim()*Npoints,normals[0].data());
00238
00239
00240 if (me == VizServer) {
00241 for (n=0; n<Npoints; n++) {
00242
00243 point_type tt;
00244 tt(0) = (sig[ n]*normals[n](0)+sig[ Npoints+n]*normals[n](1)+sig[3*Npoints+n]*normals[n](2));
00245 tt(1) = (sig[ Npoints+n]*normals[n](0)+sig[2*Npoints+n]*normals[n](1)+sig[4*Npoints+n]*normals[n](2));
00246 tt(2) = (sig[3*Npoints+n]*normals[n](0)+sig[4*Npoints+n]*normals[n](1)+sig[5*Npoints+n]*normals[n](2));
00247 st -= tt*ds[n];
00248
00249
00250 sp += p[n]*normals[n]*ds[n];
00251
00252
00253 point_type tv;
00254 tv(0) = (dev[ n]*normals[n](0)+dev[ Npoints+n]*normals[n](1)+dev[3*Npoints+n]*normals[n](2));
00255 tv(1) = (dev[ Npoints+n]*normals[n](0)+dev[2*Npoints+n]*normals[n](1)+dev[4*Npoints+n]*normals[n](2));
00256 tv(2) = (dev[3*Npoints+n]*normals[n](0)+dev[4*Npoints+n]*normals[n](1)+dev[5*Npoints+n]*normals[n](2));
00257 sv -= tv*ds[n];
00258
00259
00260 tt = -(tt-(tt(0)*normals[n](0)+tt(1)*normals[n](1)+tt(2)*normals[n](2))*normals[n]);
00261 taut[n] = tt.abs();
00262
00263 tv = -(tv-(tv(0)*normals[n](0)+tv(1)*normals[n](1)+tv(2)*normals[n](2))*normals[n]);
00264 tauv[n] = tv.abs();
00265 }
00266 }
00267 delete [] normals;
00268 delete [] ds;
00269 }
00270
00271 virtual void Output() {
00272 base::Output();
00273
00274 point_type st=0., sp=0., sv=0.;
00275 EvalForce(st,sp,sv);
00276
00277 if (MY_PROC == VizServer && SurfaceFile.c_str()[0] != '-') {
00278 char ioname[DAGHBktGFNameWidth+256];
00279 int Time = CurrentTime(base::GH(),0);
00280 std::ofstream outfile;
00281 std::ostream* out;
00282
00283 std::sprintf(ioname,"%s_%d.txt",SurfaceFile.c_str(),Time);
00284 outfile.open(ioname, std::ios::out);
00285 out = new std::ostream(outfile.rdbuf());
00286 register int n = 0;
00287 for (register int l=0; l<4*Np; l++)
00288 for (register int m=0; m<2*Np; m++) {
00289 if (m==Np) {
00290 double alpha=std::atan2(xc[n](1)-xs[1],xc[n](0)-xs[0]);
00291 if (alpha<0.) alpha+=2*pi;
00292 *out << alpha*180./pi << " " << xc[n](0) << " " << xc[n](1) << " " << xc[n](2) << " " << p[n]-p0
00293 << " " << taut[n] << " " << tauv[n] << std::endl;
00294 }
00295 n++;
00296 }
00297 outfile.close();
00298 delete out;
00299
00300 DataType fac=r0*u0*u0/2., Re=u0*2.*rd/nu;
00301 DataType fac2=2.*fac/std::sqrt(Re);
00302 std::sprintf(ioname,"%s_n_%d.txt",SurfaceFile.c_str(),Time);
00303 outfile.open(ioname, std::ios::out);
00304 out = new std::ostream(outfile.rdbuf());
00305 n = 0;
00306 for (register int l=0; l<4*Np; l++)
00307 for (register int m=0; m<2*Np; m++) {
00308 if (m==Np) {
00309 double alpha=std::atan2(xc[n](1)-xs[1],xc[n](0)-xs[0]);
00310 if (alpha<0.) alpha+=2*pi;
00311 *out << alpha*180./pi << " " << xc[n](0) << " " << xc[n](1) << " " << xc[n](2) << " "
00312 << (p[n]-p0)/fac << " " << (pa[n]-p0)/fac << " " << taut[n]/fac2 << " "
00313 << tauta[n]/fac2 << " " << tauv[n]/fac2 << " " << tauva[n]/fac2 << std::endl;
00314 }
00315 n++;
00316 }
00317 outfile.close();
00318 delete out;
00319 }
00320 }
00321
00322 virtual void Advance(double& t, double& dt) {
00323 base::Advance(t,dt);
00324
00325 if (StepsTaken(base::GH(),0)%IntegrateEvery == 0) {
00326 point_type st=0., sp=0., sv=0.;
00327 EvalForce(st,sp,sv);
00328
00329 if (MY_PROC == VizServer) {
00330 std::ofstream outfile;
00331 std::ostream* out;
00332 std::string fname = "drag.txt";
00333 outfile.open(fname.c_str(), std::ios::out | std::ios::app);
00334 out = new std::ostream(outfile.rdbuf());
00335 *out << base::t[0] << " " << st(0) << " " << st(1) << " " << st(2) << " "
00336 << sp(0) << " " << sp(1) << " " << sp(2) << " " << sv(0) << " " << sv(1) << " " << sv(2)
00337 << " " << sp(0)+sv(0) << " " << sp(1)+sv(1) << " " << sp(2)+sv(2) << std::endl;
00338 outfile.close();
00339 delete out;
00340
00341 DataType fac=r0*u0*u0*0.5*rd*rd*pi;
00342 fname = "drag_n.txt";
00343 outfile.open(fname.c_str(), std::ios::out | std::ios::app);
00344 out = new std::ostream(outfile.rdbuf());
00345 *out << base::t[0] << " " << st(0)/fac << " " << std::sqrt(st(1)*st(1)+st(2)*st(2))/fac << " "
00346 << sp(0)/fac << " " << std::sqrt(sp(1)*sp(1)+sp(2)*sp(2))/fac << " "
00347 << sv(0)/fac << " " << std::sqrt(sv(1)*sv(1)+sv(2)*sv(2))/fac << " "
00348 << (sp(0)+sv(0))/fac << " " << std::sqrt((sp(1)+sv(1))*(sp(1)+sv(1))+(sp(2)+sv(2))*(sp(2)+sv(2)))/fac << std::endl;
00349 outfile.close();
00350 delete out;
00351
00352 for (register int n=0; n<8*Np*Np; n++) {
00353 pa[n] = (pa[n]*pac+p[n])/(pac+1);
00354 tauta[n] = (tauta[n]*pac+taut[n])/(pac+1);
00355 tauva[n] = (tauva[n]*pac+tauv[n])/(pac+1);
00356 }
00357 pac++;
00358
00359 }
00360 }
00361 }
00362
00363 protected:
00364 int IntegrateEvery, Np;
00365 interpolation_type* _Interpolation;
00366 std::string SurfaceFile;
00367 double pi;
00368 DataType *dev, *sig, *p, *pa, *taut, *tauta, *tauv, *tauva;
00369 point_type *xc;
00370 int pac;
00371 };
00372
00373 #endif