File List
amroc/lbm/applications/Navier-Stokes/2d/CylinderDrag/src/Problem.h

00001 // -*- C++ -*-
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 xs[2], 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);
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 j=0; j<Ny; j++) {
00043       DataType yd = (DataType(j)+static_cast<DataType>(0.5))*dx(1) + lb(1) - xs[1];
00044       for (register int i=0; i<Nx; i++) {
00045         DataType xd = (DataType(i)+static_cast<DataType>(0.5))*dx(0) + lb(0) - xs[0];
00046         phi[i+Nx*j] = std::sqrt(xd*xd+yd*yd) - rd;
00047       }
00048     }
00049   }
00050 };
00051 
00052 
00053 class GhostFluidMethodSpecific : public GhostFluidMethod<MicroType,DIM> {
00054   typedef GhostFluidMethod<MicroType,DIM> base;
00055 public:
00056   GhostFluidMethodSpecific(boundary_type* boundary, levelset_type* levelset) : 
00057     base(boundary,levelset) {
00058     xs[0]=0.; xs[1]=0.; rd=0.; 
00059   }
00060 
00061   virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00062     base::register_at(Ctrl,prefix);
00063     RegisterAt(base::LocCtrl,"Center(1)",xs[0]);
00064     RegisterAt(base::LocCtrl,"Center(2)",xs[1]);
00065     RegisterAt(base::LocCtrl,"Radius",rd);
00066   } 
00067 };
00068 
00069  
00070 class BoundaryConditionsSpecific : public BoundaryConditionsSupplemental {
00071 public:    
00072   BoundaryConditionsSpecific(LBMType &lbm) : BoundaryConditionsSupplemental(lbm) {}
00073 };
00074 
00075 
00076 class InitialConditionSpecific : public LBMInitialCondition<LBMType,DIM> {
00077   typedef LBMInitialCondition<LBMType,DIM> base;
00078   
00079 public:    
00080   InitialConditionSpecific(LBMType &lbm) : base(lbm) {}
00081 
00082   virtual void SetupData(GridHierarchy* gh, const int& ghosts) {
00083     base::SetupData(gh,ghosts);
00084     if (base::NAux>=2) {
00085       r0 = aux[0]; u0 = aux[1];
00086     }
00087   }
00088 };
00089 
00090 
00091 class SolverSpecific : 
00092   public AMRGFMSolver<MicroType,FixupType,FlagType,DIM> {
00093   typedef AMRGFMSolver<MicroType,FixupType,FlagType,DIM> base;
00094   typedef AMRGFMInterpolation<MicroType,FixupType,FlagType,DIM> interpolation_type;
00095   typedef MicroType::InternalDataType DataType;
00096   typedef SchemeGFMFileOutput<LBMType,FixupType,FlagType,DIM> output_type;
00097   typedef interpolation_type::point_type point_type;
00098 public:
00099   SolverSpecific(IntegratorSpecific& integ, 
00100                  InitialConditionSpecific& init,
00101                  base::boundary_conditions_type& bc) : base(integ, init, bc) {
00102     SetLevelTransfer(new LBMF77LevelTransfer<LBMType,DIM>(integ.Scheme(), f_prolong, f_restrict));
00103     SetFileOutput(new SchemeGFMFileOutput<LBMType,FixupType,FlagType,DIM>(*this,integ.Scheme())); 
00104     SetFixup(new FixupSpecific(integ.Scheme()));
00105     SetFlagging(new FlaggingSpecific(*this,integ.Scheme())); 
00106     AddGFM(new GhostFluidMethodSpecific(new LBMGFMBoundary<LBMType,DIM>(integ.LBM()),
00107                                         new GFMLevelSetSpecific()));
00108     _Interpolation = new interpolation_type(*this);
00109     IntegrateEvery = 1; Np=100;
00110     SurfaceFile = "-";
00111     pi = 4.0*std::atan(1.0);
00112     dev = (DataType *)0; sig = (DataType *)0;
00113     p = (DataType *)0; pa = (DataType *)0;
00114     taut = (DataType *)0; tauta = (DataType *)0;
00115     tauv = (DataType *)0; tauva = (DataType *)0;
00116     xc = (point_type *)0;
00117     pac = 0;
00118   }  
00119 
00120   ~SolverSpecific() {
00121     DeleteGFM(_GFM[0]);
00122     delete _LevelTransfer;
00123     delete _Flagging;
00124     delete _Fixup;
00125     delete _FileOutput;
00126     delete _Interpolation;
00127     if (dev) delete [] dev;
00128     if (sig) delete [] sig;
00129     if (p) delete [] p;
00130     if (pa) delete [] pa;
00131     if (taut) delete [] taut;
00132     if (tauta) delete [] tauta;
00133     if (tauv) delete [] tauv;
00134     if (tauva) delete [] tauva;
00135     if (xc) delete [] xc;
00136   }
00137 
00138   virtual void SetupData() {
00139     base::SetupData();
00140     base::NAMRTimeSteps = 1;
00141     base::AdaptBndTimeInterpolate = 0;
00142     base::Step[0].LastTime = 1.e37;
00143     base::Step[0].VariableTimeStepping = -1;
00144     base::Step[0].dtv[0] = ((LBMIntegrator<LBMType,DIM> &)Integrator_()).LBM().TimeScale();
00145     _Interpolation->SetupData(base::PGH(), base::NGhosts());
00146 
00147     dev = new DataType[12*Np];
00148     sig = new DataType[12*Np];
00149     p = new DataType[4*Np]; pa = new DataType[4*Np];
00150     for (register int i=0; i<4*Np; i++) pa[i]=0.;
00151     taut = new DataType[4*Np]; tauta = new DataType[4*Np];
00152     for (register int i=0; i<4*Np; i++) tauta[i]=0.;
00153     tauv = new DataType[4*Np]; tauva = new DataType[4*Np];
00154     for (register int i=0; i<4*Np; i++) tauva[i]=0.;
00155     xc = new point_type[4*Np];
00156     nu = ((LBMType &)((LBMIntegrator<LBMType,DIM> &)Integrator_()).LBM()).GasViscosity();
00157     p0 = ((LBMType &)((LBMIntegrator<LBMType,DIM> &)Integrator_()).LBM()).BasePressure();
00158 
00159     std::cout << "r0=" << r0 << "   u0=" << u0 << "   rd=" << rd << "   p0=" << p0 
00160               << "   Re=" << u0*2.*rd/nu << std::endl;
00161   }
00162 
00163   virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00164     base::register_at(Ctrl,prefix);
00165     RegisterAt(base::LocCtrl,"OutputDragEvery",IntegrateEvery);
00166     RegisterAt(base::LocCtrl,"Np",Np);
00167     RegisterAt(base::LocCtrl,"SurfacePlot",SurfaceFile);
00168   } 
00169   virtual void register_at(ControlDevice& Ctrl) {
00170     base::register_at(Ctrl);
00171   }
00172 
00173   void EvalForce(point_type& st, point_type& sp, point_type& sv) {
00174     int me = MY_PROC; 
00175     int Npoints = 4*Np;
00176     DataType* dist = (DataType*) 0;    
00177     point_type* normals = new point_type[Npoints];
00178     double dalpha = 0.5*pi/Np;
00179 
00180     // Calculate equidistant points on shere's surfaces
00181     register int n=0;
00182     for (n=0; n<Npoints; n++) {
00183       double alpha = dalpha*(n+0.5);
00184       xc[n](0) = xs[0]+rd*std::cos(alpha); 
00185       xc[n](1) = xs[1]+rd*std::sin(alpha);
00186     }
00187     
00188     // Calculate pressure
00189     for (register int l=0; l<=FineLevel(base::GH()); l++) {
00190       int Time = CurrentTime(base::GH(),l);
00191       ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time, l, 
00192                                               6, base::t[l]);
00193     }
00194     // Interpolate/extrapolate pressure values and assemble results on processor VizServer
00195     _Interpolation->PointsValues(base::Work(),base::t[0],Npoints,xc,p);
00196     _Interpolation->ArrayCombine(VizServer,Npoints,p);
00197 
00198     // Calculate deviatoric stress
00199     for (register int cnt=16; cnt<=18; cnt++) {
00200       for (register int l=0; l<=FineLevel(base::GH()); l++) {
00201         int Time = CurrentTime(base::GH(),l);
00202         ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time, l, 
00203                                                 cnt, base::t[l]);
00204       }
00205       // Interpolate/extrapolate pressure values and assemble results on processor VizServer
00206       _Interpolation->PointsValues(base::Work(),base::t[0],Npoints,xc,dev+(cnt-16)*Npoints);
00207     }
00208     _Interpolation->ArrayCombine(VizServer,3*Npoints,dev);
00209  
00210     // Calculate total stress
00211     for (register int cnt=23; cnt<=25; cnt++) {
00212       for (register int l=0; l<=FineLevel(base::GH()); l++) {
00213         int Time = CurrentTime(base::GH(),l);
00214         ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time, l, 
00215                                                 cnt, base::t[l]);
00216       }
00217       // Interpolate/extrapolate pressure values and assemble results on processor VizServer
00218       _Interpolation->PointsValues(base::Work(),base::t[0],Npoints,xc,sig+(cnt-23)*Npoints);
00219     }
00220     _Interpolation->ArrayCombine(VizServer,3*Npoints,sig);
00221     
00222     // Approximate normals and assemble results on processor VizServer
00223     _Interpolation->PointsGeom(base::GFM(0).Phi(),base::t[0],Npoints,xc,normals,dist);
00224     _Interpolation->ArrayCombine(VizServer,base::Dim()*Npoints,normals[0].data());
00225 
00226     // Do the integration on the surface only on processor VizServer
00227     if (me == VizServer) {
00228       double ds = 2.*pi*rd/Npoints; 
00229       for (n=0; n<Npoints; n++) {
00230         // Total force vector
00231         point_type tt;
00232         tt(0) = (sig[        n]*normals[n](0)+sig[  Npoints+n]*normals[n](1)); 
00233         tt(1) = (sig[Npoints+n]*normals[n](0)+sig[2*Npoints+n]*normals[n](1)); 
00234         st -= tt*ds;
00235 
00236         // Hydrodynamic force vector
00237         sp += p[n]*normals[n]*ds; 
00238 
00239         // Viscous force vector
00240         point_type tv;
00241         tv(0) = (dev[        n]*normals[n](0)+dev[  Npoints+n]*normals[n](1)); 
00242         tv(1) = (dev[Npoints+n]*normals[n](0)+dev[2*Npoints+n]*normals[n](1)); 
00243         sv -= tv*ds;
00244 
00245         // Wall shear stress
00246         tt = -(tt-(tt(0)*normals[n](0)+tt(1)*normals[n](1))*normals[n]);
00247         taut[n] = tt.abs();
00248 
00249         tv = -(tv-(tv(0)*normals[n](0)+tv(1)*normals[n](1))*normals[n]);
00250         tauv[n] = tv.abs();
00251 
00252       }
00253     }
00254     delete [] normals;
00255   }
00256 
00257   virtual void Output() {
00258     base::Output();
00259 
00260     point_type st=0., sp=0., sv=0.;
00261     EvalForce(st,sp,sv);
00262     // Do the integration on the surface only on processor VizServer
00263     if (MY_PROC == VizServer && SurfaceFile.c_str()[0] != '-') {
00264       char ioname[DAGHBktGFNameWidth+256];
00265       int Time = CurrentTime(base::GH(),0);
00266       std::ofstream outfile;
00267       std::ostream* out;
00268 
00269       std::sprintf(ioname,"%s_%d.txt",SurfaceFile.c_str(),Time); 
00270       outfile.open(ioname, std::ios::out);
00271       out = new std::ostream(outfile.rdbuf());
00272       for (register int n=0; n<4*Np; n++) {
00273         double alpha=std::atan2(xc[n](1)-xs[1],xc[n](0)-xs[0]);
00274         if (alpha<0.) alpha+=2*pi;        
00275         *out << alpha*180./pi << " " << xc[n](0) << " " << xc[n](1) << " " << p[n]-p0 << " " << pa[n]-p0 
00276              << " " << taut[n] << " " << tauta[n] << " " << tauv[n] << " " << tauva[n] << std::endl;
00277       }
00278       outfile.close();
00279       delete out;   
00280 
00281       DataType fac=r0*u0*u0/2., Re=u0*2.*rd/nu;
00282       DataType fac2=2.*fac/std::sqrt(Re);
00283       std::sprintf(ioname,"%s_n_%d.txt",SurfaceFile.c_str(),Time); 
00284       outfile.open(ioname, std::ios::out);
00285       out = new std::ostream(outfile.rdbuf());
00286       for (register int n=0; n<4*Np; n++) {
00287         double alpha=std::atan2(xc[n](1)-xs[1],xc[n](0)-xs[0]);
00288         if (alpha<0.) alpha+=2*pi;        
00289         *out << alpha*180./pi << " " << xc[n](0) << " " << xc[n](1) << " " << (p[n]-p0)/fac << " " << (pa[n]-p0)/fac 
00290              << " " << taut[n]/fac2 << " " << tauta[n]/fac2 << " " << tauv[n]/fac2 << " " << tauva[n]/fac2 << std::endl;
00291       }
00292       outfile.close();
00293       delete out;   
00294     }
00295   }
00296       
00297   virtual void Advance(double& t, double& dt) {
00298     base::Advance(t,dt);
00299 
00300     if (StepsTaken(base::GH(),0)%IntegrateEvery == 0) {
00301       point_type st=0., sp=0., sv=0.;
00302       EvalForce(st,sp,sv);
00303       // Append to output file only on processor VizServer
00304       if (MY_PROC == VizServer) {
00305         std::ofstream outfile;
00306         std::ostream* out;
00307         std::string fname = "drag.txt";
00308         outfile.open(fname.c_str(), std::ios::out | std::ios::app);
00309         out = new std::ostream(outfile.rdbuf());
00310         *out << base::t[0] << " " << st(0) << " " << st(1) << " " << sp(0) << " " << sp(1) 
00311              << " " << sv(0) << " " << sv(1) << " " << sp(0)+sv(0) << " " << sp(1)+sv(1) << std::endl;
00312         outfile.close();
00313         delete out;   
00314 
00315         DataType fac=r0*u0*u0*rd;
00316         fname = "drag_n.txt";
00317         outfile.open(fname.c_str(), std::ios::out | std::ios::app);
00318         out = new std::ostream(outfile.rdbuf());
00319         *out << base::t[0] << " " << st(0)/fac << " " << st(1)/fac << " " << sp(0)/fac << " " << sp(1)/fac 
00320              << " " << sv(0)/fac << " " << sv(1)/fac << " " << (sp(0)+sv(0))/fac << " " << (sp(1)+sv(1))/fac << std::endl;
00321         outfile.close();
00322         delete out;   
00323 
00324         for (register int n=0; n<4*Np; n++) {
00325           pa[n] = (pa[n]*pac+p[n])/(pac+1);
00326           tauta[n] = (tauta[n]*pac+taut[n])/(pac+1);
00327           tauva[n] = (tauva[n]*pac+tauv[n])/(pac+1);
00328         }
00329         pac++;
00330       }
00331     }
00332   }
00333 
00334 protected:
00335   int IntegrateEvery, Np;
00336   interpolation_type* _Interpolation; 
00337   std::string SurfaceFile;
00338   double pi;
00339   DataType *dev, *sig, *p, *pa, *taut, *tauta, *tauv, *tauva;
00340   point_type *xc;
00341   int pac;
00342 }; 
00343 
00344 #endif
00345 
00346 
00347
Generated on Thu Jun 30 02:20:12 2016 for AMROC Applications by 1.6.3