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