00001
00002
00003 #ifndef AMROC_FLUIDPROBLEM_H
00004 #define AMROC_FLUIDPROBLEM_H
00005 #define DIM 3
00006
00007 #define MaxIntPoints (20)
00008 #define MaxIntPoints_init (3)
00009
00010 #include "LBMProblem.h"
00011 #include "LBMD3Q19.h"
00012
00013 typedef double DataType;
00014 typedef LBMD3Q19<DataType> LBMType;
00015
00016 #define OWN_LBMSCHEME
00017 #define OWN_AMRSOLVER
00018 #define OWN_INITIALCONDITION
00019 #define OWN_BOUNDARYCONDITION
00020 #include "LBMStdProblem.h"
00021 #include "F77Interfaces/F77GFMBoundary.h"
00022 #include "AMRELCGFMSolver.h"
00023 #include "LBMGFMBoundary.h"
00024 #include "Interfaces/SchemeGFMFileOutput.h"
00025
00026
00027 #define PATH_LOG
00028 #define maxPathlines 3
00029 #define maxPathLength 1000
00030 #define maxPointClouds 3
00031 #define tol 1.0e-5
00032 #define MAXVEL 100.0
00033
00034 class PathPointSourceCtrl : public controlable {
00035 typedef ads::FixedArray<8,DataType> PVType;
00036 public:
00037 PathPointSourceCtrl() {
00038 path.resize(0);
00039 path.reserve(maxPathLength);
00040 }
00041
00042 virtual ~PathPointSourceCtrl() {
00043
00044 }
00045
00046 virtual void register_at(ControlDevice& Ctrl) {}
00047 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00048 PLocCtrl = Ctrl.getSubDevice(prefix);
00049 RegisterAt(PLocCtrl,"point(0)",pointSource[0]);
00050 RegisterAt(PLocCtrl,"point(1)",pointSource[1]);
00051 RegisterAt(PLocCtrl,"point(2)",pointSource[2]);
00052 }
00053
00054 mutable PVType pointSource, pointSource_old;
00055 mutable std::vector<PVType> path;
00056 ControlDevice PLocCtrl;
00057 };
00058
00059 class LBMSpecific : public LBMType {
00060 typedef LBMType base;
00061 typedef ads::FixedArray<8,DataType> PVType;
00062 typedef ads::FixedArray<10,DataType> PVCType;
00063 public:
00064 LBMSpecific() : base(), omega(1.), u_p_avg(1.), l0_p(-1.) {
00065
00066 }
00067
00068 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00069 base::register_at(Ctrl,prefix);
00070 RegisterAt(base::LocCtrl,"Omega",omega);
00071 RegisterAt(base::LocCtrl,"u_p_avg",u_p_avg);
00072 RegisterAt(base::LocCtrl,"l0_p",l0_p);
00073
00074 }
00075
00076 virtual void SetupData(GridHierarchy* gh, const int& ghosts) {
00077 base::SetupData(gh,ghosts);
00078
00079 const double* bndry = base::GH().wholebndry();
00080 const int nbndry = base::GH().nbndry();
00081 if (l0_p <= 0. ) {
00082 l0_p=bndry[nbndry*(1+2*1)]-bndry[nbndry*(0+2*1)];
00083 if (nbndry>1) l0_p=bndry[nbndry*(1+2*1)]-bndry[1+nbndry*(1+2*1)];
00084 }
00085 DataType omega = base::Omega(base::TimeScale());
00086
00087
00088 double Re_p=u_p_avg*l0_p/base::GasViscosity();
00089 double Re_lbm=(u_p_avg/base::VelocityScale())*(l0_p/base::LengthScale())/base::LatticeViscosity(omega);
00090 ( comm_service::log() << "l0_p=" << l0_p << " u_p_avg=" << u_p_avg
00091 << " Re_p=" << Re_p << " Re_lbm=" << Re_lbm
00092 << "\n omega=" << omega
00093 << " SpeedUp=" << base::SpeedUp()
00094 << " u_avg_LBM=" << u_p_avg/base::VelocityScale()
00095 << " SpeedRatio=" << u_p_avg/base::VelocityScale()*std::sqrt(DataType(3.))
00096 << " base::TO()=" << base::T0()
00097 << "\n Re_p-Re_lbm=" << Re_p-Re_lbm
00098 << std::endl ).flush();
00099
00100 if ( std::fabs(Re_p-Re_lbm)>DBL_EPSILON*LB_FACTOR*LB_FACTOR*LB_FACTOR ) {
00101 std::cerr << "Physical scaling requires fabs(Re_p-Re_lbm)<DBL_EPSILON*LB_FACTOR*LB_FACTOR*LB_FACTOR !" << std::endl;
00102 std::cerr << "Physical scaling requires fabs(Re_p-Re_lbm)<DBL_EPSILON*LB_FACTOR*LB_FACTOR*LB_FACTOR !" << std::endl;
00103 std::cerr << "Re_p " << Re_p << " Re_lbm " << Re_lbm << " std::fabs(Re_p-Re_lbm) " << std::fabs(Re_p-Re_lbm) << std::endl;
00104 std::exit(-1);
00105 }
00106 if (omega<0.5 || omega>2.) {
00107 std::cerr << "0.5 <= omega <= 2 required." << std::endl;
00108 std::exit(-1);
00109 }
00110 base::SetGas(base::GasDensity(),base::GasViscosity(),base::GasSpeedofSound());
00111 ( comm_service::log() << "FluidProblem SetupData() D3Q19: Gas_rho=" << base::GasDensity() << " Gas_Cs=" << base::GasSpeedofSound()
00112 << " Gas_nu=" << base::GasViscosity() << std::endl ).flush();
00113
00114 }
00115
00116 void updateCloud(int n, PVType val) const {}
00117
00118 void updatePath(int n, PVType val) const {
00119 PointSourceCtrl[n].path.push_back(val);
00120 #ifdef PATH_DEBUG
00121 ( comm_service::log() << " PointSourceCtrl[" << n << "]." << " path[" << PointSourceCtrl[n].path.size()-1 << "]="
00122 << PointSourceCtrl[n].path[PointSourceCtrl[n].path.size()-1] << std::endl ).flush();
00123 #endif
00124 }
00125
00126 void ouputPath(int n, DataType t) const {
00127 char fname[256];
00128 sprintf(fname,"pathLine_%02d_%05d.vtk",n,StepsTaken((GridHierarchy&) base::GH(),0) );
00129 std::ofstream ofs;
00130 ofs.open(fname, std::ios::out);
00131
00132 ofs << "# vtk DataFile Version 2.0\n" << "pathLine Test"
00133 << "\nASCII\nDATASET UNSTRUCTURED_GRID\nFIELD FieldData " << 1 << std::endl;
00134 ofs << "TIME 1 1 float\n " << t << "\nPOINTS " << PointSourceCtrl[n].path.size() << " float\n";
00135 for (unsigned int i = 0; i < PointSourceCtrl[n].path.size(); i++ ) {
00136 ofs << PointSourceCtrl[n].path[i][0] << " " << PointSourceCtrl[n].path[i][1] << " " << PointSourceCtrl[n].path[i][2] << std::endl;
00137
00138 }
00139 ofs << "CELLS " << PointSourceCtrl[n].path.size()-1 << " " << (PointSourceCtrl[n].path.size()-1)*3 << std::endl;
00140 for (unsigned int i = 0; i < PointSourceCtrl[n].path.size()-1; i++ ) {
00141 ofs << 2 << " " << i << " " << i+1 << std::endl;
00142 }
00143 ofs << "CELL_TYPES " << PointSourceCtrl[n].path.size()-1 << std::endl;
00144 for (unsigned int i = 0; i < PointSourceCtrl[n].path.size()-1; i++ )
00145 ofs << 3 << std::endl;
00146 ofs << "CELL_DATA " << PointSourceCtrl[n].path.size()-1 << "\nFIELD Data " << 5 << std::endl;
00147 ofs << "rho " << 1 << " " << PointSourceCtrl[n].path.size()-1 << " double\n";
00148 for (unsigned int i = 1; i < PointSourceCtrl[n].path.size(); i++ ) {
00149 ofs << DataType(0.5)*(PointSourceCtrl[n].path[i-1][3]+PointSourceCtrl[n].path[i][3]) << std::endl;
00150 }
00151 ofs << "u " << 1 << " " << PointSourceCtrl[n].path.size()-1 << " double\n";
00152 for (unsigned int i = 1; i < PointSourceCtrl[n].path.size(); i++ ) {
00153 ofs << DataType(0.5)*(PointSourceCtrl[n].path[i-1][4]+PointSourceCtrl[n].path[i][4]) << std::endl;
00154 }
00155 ofs << "v " << 1 << " " << PointSourceCtrl[n].path.size()-1 << " double\n";
00156 for (unsigned int i = 1; i < PointSourceCtrl[n].path.size(); i++ ) {
00157 ofs << DataType(0.5)*(PointSourceCtrl[n].path[i-1][5]+PointSourceCtrl[n].path[i][5]) << std::endl;
00158 }
00159 ofs << "w " << 1 << " " << PointSourceCtrl[n].path.size()-1 << " double\n";
00160 for (unsigned int i = 1; i < PointSourceCtrl[n].path.size(); i++ ) {
00161 ofs << DataType(0.5)*(PointSourceCtrl[n].path[i-1][6]+PointSourceCtrl[n].path[i][6]) << std::endl;
00162 }
00163 ofs << "P " << 1 << " " << PointSourceCtrl[n].path.size()-1 << " double\n";
00164 for (unsigned int i = 1; i < PointSourceCtrl[n].path.size(); i++ ) {
00165 ofs << DataType(0.5)*(PointSourceCtrl[n].path[i-1][7]+PointSourceCtrl[n].path[i][7]) << std::endl;
00166 }
00167 ofs.close();
00168 }
00169
00170 private:
00171 DataType omega, u_p_avg, l0_p;
00172
00173 protected:
00174 int track_every, num_pathLines, num_clouds;
00175 PathPointSourceCtrl PointSourceCtrl[maxPathlines];
00176 };
00177
00178 class InitialConditionSpecific : public LBMInitialCondition<LBMType,DIM> {
00179 typedef LBMInitialCondition<LBMType,DIM> base;
00180 public:
00181 InitialConditionSpecific(LBMType &lbm) : base(lbm) {}
00182 typedef base::MicroType MicroType;
00183 typedef base::MacroType MacroType;
00184
00185 virtual void SetGrid(vec_grid_data_type& fvec, grid_data_type& workvec, const int& level) {
00186 MacroType q, q_;
00187 DataType elevation = 15., profile, height, ldx;
00188 BBox wb = base::GH().wholebbox();
00189 BBox bb = fvec.bbox();
00190 DCoords wc;
00191
00192 q_(0) = aux[0]/LBM().DensityScale();
00193 q_(1) = aux[1]/LBM().VelocityScale();
00194 q_(2) = aux[2]/LBM().VelocityScale();
00195 q_(3) = aux[3]/LBM().VelocityScale();
00196
00197 MicroType *f = (MicroType *)fvec.databuffer();
00198
00199 int Nx = fvec.extents(0), Ny = fvec.extents(1);
00200 int b = fvec.bottom();
00201 assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1));
00202 int mxs = std::max(bb.lower(0)/bb.stepsize(0),fvec.lower(0)/fvec.stepsize(0));
00203 int mys = std::max(bb.lower(1)/bb.stepsize(1),fvec.lower(1)/fvec.stepsize(1));
00204 int mzs = std::max(bb.lower(2)/bb.stepsize(2),fvec.lower(2)/fvec.stepsize(2));
00205 int mxe = std::min(bb.upper(0)/bb.stepsize(0),fvec.upper(0)/fvec.stepsize(0));
00206 int mye = std::min(bb.upper(1)/bb.stepsize(1),fvec.upper(1)/fvec.stepsize(1));
00207 int mze = std::min(bb.upper(2)/bb.stepsize(2),fvec.upper(2)/fvec.stepsize(2));
00208 int lc_indx[3];
00209
00210 lc_indx[0] = base::GH().worldCoords(fvec.lower(),fvec.stepsize())(0);
00211 lc_indx[1] = base::GH().worldCoords(fvec.lower(),fvec.stepsize())(1);
00212
00213 ldx = DataType(fvec.stepsize(2))/DataType(wb.stepsize(2))*LBM().L0();
00214
00215 for (register int k=mzs; k<=mze; k++) {
00216 lc_indx[2] = (int) std::floor(DataType(k)*DataType(ldx));
00217 wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
00218 height = DataType(wc(2))*DataType(wb.stepsize(2));
00219
00220
00221 if ( height < elevation ) {
00222 profile = DataType(height*height)/DataType(elevation*elevation);
00223
00224 q(0) = q_(0);
00225 q(1) = q_(1)*profile;
00226 q(2) = q_(2)*profile;
00227 q(3) = q_(3)*profile;
00228 }
00229 else {
00230 q = q_;
00231
00232 }
00233 for (register int j=mys; j<=mye; j++)
00234 for (register int i=mxs; i<=mxe; i++)
00235 f[b+LBM().idx(i,j,k,Nx,Ny)] = LBM().Equilibrium(q);
00236 }
00237
00238 }
00239 };
00240
00241 class BoundaryConditionsSpecific : public LBMBoundaryConditions<LBMType,DIM> {
00242 typedef LBMBoundaryConditions<LBMType,DIM> base;
00243 public:
00244 typedef base::MicroType MicroType;
00245 typedef base::MacroType MacroType;
00246 BoundaryConditionsSpecific(LBMType &lbm) : base(lbm) {}
00247
00248 DataType caux[6];
00249 int ctype, cside, cscaling, ncaux;
00250
00251 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00252 base::register_at(Ctrl, prefix);
00253 base::LocCtrl = Ctrl.getSubDevice(prefix+"BoundaryConditionCustom");
00254 RegisterAt(base::LocCtrl,"Type",ctype);
00255 RegisterAt(base::LocCtrl,"Side",cside);
00256 RegisterAt(base::LocCtrl,"Scaling",cscaling);
00257 RegisterAt(base::LocCtrl,"NCAux",ncaux);
00258 char VariableName[16];
00259 for (int i=0; i<6; i++) {
00260 caux[i]=DataType(1.);
00261 std::sprintf(VariableName,"CAux(%d)",i+1);
00262 RegisterAt(base::LocCtrl,VariableName,caux[i]);
00263 }
00264 }
00265
00266 virtual void SetupData(GridHierarchy* gh, const int& ghosts) {
00267 base::SetupData(gh,ghosts);
00268 const double* bndry = base::GH().wholebndry();
00269 const int nbndry = base::GH().nbndry();
00270 DataType l0_p=bndry[nbndry*(1+2*1)]-bndry[nbndry*(0+2*1)];
00271 if (nbndry>1) l0_p=bndry[nbndry*(1+2*1)]-bndry[1+nbndry*(1+2*1)];
00272 DataType omega = LBM().Omega(LBM().TimeScale());
00273
00274 DataType u_p_avg=caux[1];
00275 double Re_p=u_p_avg*l0_p/LBM().GasViscosity();
00276 double Re_lbm=(u_p_avg/LBM().VelocityScale())*(l0_p/LBM().LengthScale())/LBM().LatticeViscosity(omega);
00277 ( comm_service::log() << "l0_p=" << l0_p << " u_p_avg=" << u_p_avg
00278 << " Re_p=" << Re_p << " Re_lbm=" << Re_lbm
00279 << "\n omega=" << omega
00280 << " SpeedUp=" << LBM().SpeedUp()
00281 << " u_avg_LBM=" << u_p_avg/LBM().VelocityScale()
00282 << " SpeedRatio=" << u_p_avg/LBM().VelocityScale()*std::sqrt(DataType(3.))
00283 << " LBM().TO()=" << LBM().T0()
00284 << "\n Re_p-Re_lbm=" << Re_p-Re_lbm
00285 << std::endl ).flush();
00286 assert(std::fabs(Re_p-Re_lbm)<DBL_EPSILON*LB_FACTOR*LB_FACTOR*LB_FACTOR);
00287 }
00288
00289 virtual void SetBndry(vec_grid_data_type &fvec, const int& level, const BBox &bb,
00290 const int &dir, const double& time) {
00291
00292
00293
00294 base::SetBndry(fvec,level,bb,dir,time);
00295 BBox wholebbox = base::GH().wholebbox();
00296 BBox inside = bb*coarsen(wholebbox,wholebbox.stepsize()/wholebbox.stepsize());
00297 if (!inside.empty())
00298 LBM().BCStandard(fvec,bb,LBMType::NoSlipWall,dir);
00299
00300 int scaling = LBMType::Physical;
00301 int Nx = fvec.extents(0), Ny = fvec.extents(1);
00302 int b = fvec.bottom();
00303 assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1));
00304 int mys = std::max(bb.lower(1)/bb.stepsize(1),fvec.lower(1)/fvec.stepsize(1));
00305 int mzs = std::max(bb.lower(2)/bb.stepsize(2),fvec.lower(2)/fvec.stepsize(2));
00306 int mxe = std::min(bb.upper(0)/bb.stepsize(0),fvec.upper(0)/fvec.stepsize(0));
00307 int mye = std::min(bb.upper(1)/bb.stepsize(1),fvec.upper(1)/fvec.stepsize(1));
00308 int mze = std::min(bb.upper(2)/bb.stepsize(2),fvec.upper(2)/fvec.stepsize(2));
00309 MicroType *f = (MicroType *)fvec.databuffer();
00310
00311 DataType rho=caux[0], a0=caux[1], a1=caux[2], a2=caux[3], elevation=caux[4], jr, hr;
00312
00313 MacroType q, q_;
00314 DataType profile, height, ldx;
00315
00316 if (scaling==LBMType::Physical) {
00317 rho /= LBM().DensityScale();
00318 a0 /= LBM().VelocityScale();
00319 a1 /= LBM().VelocityScale();
00320 a2 /= LBM().VelocityScale();
00321 elevation /= LBM().L0();
00322 }
00323
00324 if (bb.upper(0)+bb.stepsize(0)==wholebbox.lower(0)) {
00325
00326
00327 BBox wb = base::GH().wholebbox();
00328 DCoords wc;
00329 q_(0) = rho;
00330 q_(1) = a0;
00331 q_(2) = a1;
00332 q_(3) = a2;
00333
00334 int lc_indx[3];
00335
00336 lc_indx[0] = base::GH().worldCoords(fvec.lower(),fvec.stepsize())(0);
00337 lc_indx[1] = base::GH().worldCoords(fvec.lower(),fvec.stepsize())(1);
00338 ldx = DataType(fvec.stepsize(2))/DataType(wb.stepsize(2))*DataType(LBM().L0());
00339
00340 for (register int k=mzs; k<=mze; k++) {
00341 lc_indx[2] = k;
00342
00343 wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
00344 height = DataType(wc(2))*DataType(wb.stepsize(2));
00345 if ( height < elevation ) {
00346 profile = DataType(height*height)/DataType(elevation*elevation);
00347
00348 q(0) = q_(0);
00349 q(1) = q_(1)*profile;
00350 q(2) = q_(2)*profile;
00351 q(3) = q_(3)*profile;
00352 }
00353 else {
00354 q = q_;
00355
00356 }
00357 for (register int j=mys; j<=mye; j++) {
00358 f[b+LBM().idx(mxe,j,k,Nx,Ny)] = LBM().Equilibrium(q);
00359 }
00360
00361
00362 }
00363 }
00364 }
00365 };
00366
00367
00368 template <class LBMType, int dim>
00369 class LBMELCGFMBoundary : public LBMGFMBoundary<LBMType,dim> {
00370 typedef typename LBMType::MicroType MicroType;
00371 typedef typename MicroType::InternalDataType DataType;
00372 typedef LBMGFMBoundary<LBMType,dim> base;
00373 public:
00374 typedef typename base::vec_grid_data_type vec_grid_data_type;
00375 typedef typename base::grid_data_type grid_data_type;
00376 typedef typename base::point_type point_type;
00377 typedef AMRELCGFMSolver<MicroType,FixupType,FlagType,dim> amr_elc_solver;
00378
00379 LBMELCGFMBoundary(LBMType &scheme, amr_elc_solver& solver) : base(scheme), _solver(solver) {
00380 base::SetNAux(base::Dim());
00381 }
00382
00383 virtual void SetBndryAux(vec_grid_data_type& gdu, grid_data_type& gdphi,
00384 const MicroType* u, DataType* aux, const int& Level,
00385 double t, const int& nc, const int* idx,
00386 const point_type* xc, const DataType* distance,
00387 const point_type* normal) {
00388 _solver.ExtractBoundaryVelocities(gdu.bbox(),nc,xc,reinterpret_cast<point_type *>(aux));
00389 }
00390
00391 protected:
00392 amr_elc_solver& _solver;
00393 };
00394
00395 class SolverObjects {
00396 public:
00397 SolverObjects() {
00398 _LBMSpecific = new LBMSpecific();
00399 _IntegratorSpecific = new IntegratorSpecific(*_LBMSpecific);
00400 _InitialConditionSpecific = new InitialConditionSpecific(*_LBMSpecific);
00401 _BoundaryConditionsSpecific = new BoundaryConditionsSpecific(*_LBMSpecific);
00402 }
00403
00404 ~SolverObjects() {
00405 delete _BoundaryConditionsSpecific;
00406 delete _InitialConditionSpecific;
00407 delete _IntegratorSpecific;
00408 delete _LBMSpecific;
00409 }
00410
00411 protected:
00412 LBMSpecific *_LBMSpecific;
00413 IntegratorSpecific *_IntegratorSpecific;
00414 InitialConditionSpecific *_InitialConditionSpecific;
00415 BoundaryConditionsSpecific *_BoundaryConditionsSpecific;
00416 };
00417
00418 class PointCloudCtrl2 : public controlable {
00419 typedef ads::FixedArray<3,DataType> PType;
00420 typedef ads::FixedArray<10,DataType> PVCType;
00421 typedef AMRInterpolation<MicroType,DIM> interpolation_type;
00422 typedef interpolation_type::point_type point_type;
00423 public:
00424 PointCloudCtrl2() {
00425 emit_time = 0.;
00426 _IntName = "ptrack.txt";
00427 num_points = 0;
00428 pointCloud.resize(0);
00429 }
00430
00431 virtual ~PointCloudCtrl2() {}
00432
00433 virtual void register_at(ControlDevice& Ctrl) {}
00434 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00435 CLocCtrl = Ctrl.getSubDevice(prefix);
00436 char VariableName[32];
00437 for (int nc=0; nc<MaxIntPoints_init; nc++) {
00438 for (int d=0; d<DIM; d++) {
00439 std::sprintf(VariableName,"Point(%d,%d)",nc+1,d+1);
00440 RegisterAt(CLocCtrl,VariableName,_IntPoints_init[nc](d));
00441 }
00442 }
00443 RegisterAt(CLocCtrl,"FileName",_IntName);
00444 RegisterAt(CLocCtrl,"type",type);
00445 RegisterAt(CLocCtrl,"num_initPoints",num_initPoints);
00446 RegisterAt(CLocCtrl,"emit_init",emit_init);
00447 RegisterAt(CLocCtrl,"emit_intrvl",emit_intrvl);
00448 }
00449
00450 virtual void Restart(std::ifstream& ifs, int& pos) {
00451 ifs.seekg(pos);
00452 num_points = 0;
00453 emit_count = -1;
00454 ifs.read((char*)&num_points,sizeof(int));
00455 ifs.read((char*)&emit_count,sizeof(int));
00456 ifs.read((char*)&emit_time,sizeof(DataType));
00457 ifs.read((char*)&_IntPoints_init,num_initPoints*sizeof(point_type));
00458 ifs.read((char*)&_IntPoints,num_points*sizeof(point_type));
00459 ifs.read((char*)&age,num_points*sizeof(DataType));
00460
00461 for (int k=0; k<num_points; k++) {
00462 ( comm_service::log() << k << " " << _IntPoints[k] << std::endl ).flush();
00463 }
00464
00465 }
00466
00467 virtual void Checkpointing(std::ofstream& ofs) {
00468 ( comm_service::log() << "\nCHECKPOINTING CloudCtrl " << std::endl ).flush();
00469 ( comm_service::log() << "num_points " << num_points << std::endl ).flush();
00470 ( comm_service::log() << "emit_count " << emit_count << std::endl ).flush();
00471 ofs.write((char*)&num_points,sizeof(int));
00472 ofs.write((char*)&emit_count,sizeof(int));
00473 ofs.write((char*)&emit_time,sizeof(DataType));
00474 ofs.write((char*)&_IntPoints_init,num_initPoints*sizeof(point_type));
00475 ofs.write((char*)&_IntPoints,num_points*sizeof(point_type));
00476 ofs.write((char*)&age,num_points*sizeof(DataType));
00477 }
00478
00479 int type, num_initPoints;
00480 mutable int cloudUpdate, emit_count, num_points;
00481 DataType emit_init, emit_intrvl;
00482 mutable DataType emit_time;
00483
00484 mutable std::vector<PVCType> pointCloud, pointCloud_old;
00485 mutable std::vector<int> cloudCheck;
00486 ControlDevice CLocCtrl;
00487
00488 mutable point_type _IntPoints_init[MaxIntPoints_init];
00489 mutable point_type _IntPoints[MaxIntPoints];
00490 mutable DataType age[MaxIntPoints];
00491 std::string _IntName;
00492 };
00493
00494 class FluidSolverSpecific : protected SolverObjects,
00495 public AMRELCGFMSolver<MicroType,FixupType,FlagType,DIM> {
00496 typedef MicroType::InternalDataType DataType;
00497 typedef AMRELCGFMSolver<MicroType,FixupType,FlagType,DIM> base;
00498 typedef AMRInterpolation<MicroType,DIM> interpolation_type;
00499 typedef interpolation_type::point_type point_type;
00500 typedef SchemeGFMFileOutput<LBMType,FixupType,FlagType,DIM> output_type;
00501 typedef ads::FixedArray<10,DataType> PVCType;
00502 public:
00503 FluidSolverSpecific() : SolverObjects(), base(*SolverObjects::_IntegratorSpecific,
00504 *SolverObjects::_InitialConditionSpecific,
00505 *SolverObjects::_BoundaryConditionsSpecific) {
00506 SetLevelTransfer(new LBMF77LevelTransfer<LBMType,DIM>(SolverObjects::_IntegratorSpecific->Scheme(), f_prolong, f_restrict));
00507 SetFileOutput(new output_type(*this,SolverObjects::_IntegratorSpecific->Scheme()));
00508
00509 SetFixup(new FixupSpecific(SolverObjects::_IntegratorSpecific->Scheme()));
00510 SetFlagging(new FlaggingSpecific(*this,SolverObjects::_IntegratorSpecific->Scheme()));
00511 AddGFM(new GhostFluidMethod<MicroType,DIM>(new LBMELCGFMBoundary<LBMType,DIM>(SolverObjects::_IntegratorSpecific->LBM(),*this),
00512 new CPTLevelSet<DataType,DIM>()));
00513
00514 SetCoupleGFM(0);
00515
00516 _Interpolation = new interpolation_type();
00517 }
00518
00519 ~FluidSolverSpecific() {
00520 DeleteGFM(_GFM[0]);
00521 delete _LevelTransfer;
00522 delete _Flagging;
00523 delete _Fixup;
00524 delete _FileOutput;
00525 delete _BoundaryConditionsSpecific;
00526 _BoundaryConditionsSpecific = (BoundaryConditionsSpecific *) 0;
00527 delete _InitialConditionSpecific;
00528 _InitialConditionSpecific = (InitialConditionSpecific *) 0;
00529 delete _IntegratorSpecific;
00530 _IntegratorSpecific = (IntegratorSpecific *) 0;
00531 delete _Interpolation;
00532 _Interpolation = (AMRELCGFMSolver<MicroType,FixupType,FlagType,DIM>::interpolation_type *) 0;
00533 }
00534
00535 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00536 base::register_at(Ctrl,prefix);
00537 RegisterAt(base::LocCtrl,"num_clouds",num_clouds);
00538 RegisterAt(base::LocCtrl,"track_every",track_every);
00539 for (register int j=0; j < maxPointClouds; j++) {
00540 char VariableName1[32];
00541 std::sprintf(VariableName1,"PointCloud(%02d)",j+1);
00542 CloudCtrl[j].register_at(base::LocCtrl,std::string(VariableName1));
00543 }
00544
00545 }
00546 virtual void register_at(ControlDevice& Ctrl) {
00547 base::register_at(Ctrl);
00548 }
00549
00550 virtual void SendBoundaryData() {
00551 START_WATCH
00552 for (register int l=0; l<=FineLevel(base::GH()); l++) {
00553 int Time = CurrentTime(base::GH(),l);
00554 ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time,
00555 l, base::Dim()+4, base::t[l]);
00556 if (CurrentTime(base::GH(),base::CouplingLevel) != Time)
00557 ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time+TimeStep(base::U(),l),
00558 l, base::Dim()+4, base::t[l]+base::dt[l]);
00559 }
00560 END_WATCH(FLUID_CPL_PRESSURE_CALCULATE)
00561 base::SendBoundaryData();
00562 }
00563
00564 virtual void SetupData() {
00565 ( comm_service::log() << " fluid setup begun" << std::endl ).flush();
00566 base::SetupData();
00567 base::AdaptBndTimeInterpolate = 0;
00568 base::NAMRTimeSteps = 1;
00569 base::Step[0].LastTime = 1.e37;
00570 base::Step[0].VariableTimeStepping = -1;
00571 base::Step[0].dtv[0] = ((LBMIntegrator<LBMType,DIM> &)Integrator_()).LBM().TimeScale();
00572
00573 _Interpolation->SetupData(base::PGH(), base::NGhosts());
00574
00575 for (register int c = 0; c < num_clouds; c++) {
00576 CloudCtrl[c].num_points = CloudCtrl[c].num_initPoints;
00577 for (int m=0; m<CloudCtrl[c].num_points;m++) {
00578 CloudCtrl[c]._IntPoints[m] = CloudCtrl[c]._IntPoints_init[m];
00579 CloudCtrl[c].age[m] = 0;
00580 }
00581 }
00582
00583 }
00584
00585 virtual void Advance(double& t, double& dt) {
00586 base::Advance(t,dt);
00587
00588 int me = MY_PROC;
00589
00590 int nsteps = StepsTaken((GridHierarchy&) base::GH(),0);
00591
00592
00593 for (register int c = 0; c < num_clouds; c++) {
00594 if (t < CloudCtrl[c].emit_init ) return;
00595 if (CloudCtrl[c].num_points<=0 || CloudCtrl[c].num_points>MaxIntPoints) break;
00596 if (CloudCtrl[c].num_points + CloudCtrl[c].num_initPoints < MaxIntPoints ) {
00597 if ( (modulus((t - CloudCtrl[c].emit_init), CloudCtrl[c].emit_intrvl) < dt) && nsteps > 1 ) {
00598 if ( (t > (CloudCtrl[c].emit_time + dt) ) ) {
00599 for (int p = 0; p < CloudCtrl[c].num_initPoints; p++ ) {
00600 CloudCtrl[c]._IntPoints[CloudCtrl[c].num_points](0) = CloudCtrl[c]._IntPoints_init[p](0);
00601 CloudCtrl[c]._IntPoints[CloudCtrl[c].num_points](1) = CloudCtrl[c]._IntPoints_init[p](1);
00602 CloudCtrl[c]._IntPoints[CloudCtrl[c].num_points](2) = CloudCtrl[c]._IntPoints_init[p](2);
00603 CloudCtrl[c].age[CloudCtrl[c].num_points] = 0;
00604 CloudCtrl[c].num_points++;
00605
00606 }
00607 CloudCtrl[c].emit_time = t;
00608 CloudCtrl[c].emit_count++;
00609 }
00610 }
00611 }
00612 }
00613
00614 for (register int c = 0; c < num_clouds; c++) {
00615 const int numP = CloudCtrl[c].num_points;
00616 DataType* rho = new DataType[numP];
00617 DataType* u = new DataType[numP];
00618 DataType* v = new DataType[numP];
00619 DataType* w = new DataType[numP];
00620 DataType* p = new DataType[numP];
00621
00622 int Time, press_cnt;
00623
00624
00625 for (register int l=0; l<=FineLevel(base::GH()); l++) {
00626 Time = CurrentTime(base::GH(),l);
00627 press_cnt = 7;
00628 ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time, l,
00629 press_cnt, base::t[l]);
00630 }
00631 _Interpolation->PointsValuesPar(base::Work(),base::t[0],numP,CloudCtrl[c]._IntPoints,p);
00632 _Interpolation->ArrayCombine(VizServer,numP,p);
00633
00634 for (register int l=0; l<=FineLevel(base::GH()); l++) {
00635 Time = CurrentTime(base::GH(),l);
00636 press_cnt = 2;
00637 ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time, l,
00638 press_cnt, base::t[l]);
00639 }
00640 _Interpolation->PointsValuesPar(base::Work(),base::t[0],numP,CloudCtrl[c]._IntPoints,u);
00641 _Interpolation->ArrayCombine(VizServer,numP,u);
00642
00643 for (register int l=0; l<=FineLevel(base::GH()); l++) {
00644 Time = CurrentTime(base::GH(),l);
00645 press_cnt = 3;
00646 ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time, l,
00647 press_cnt, base::t[l]);
00648 }
00649 _Interpolation->PointsValuesPar(base::Work(),base::t[0],numP,CloudCtrl[c]._IntPoints,v);
00650 _Interpolation->ArrayCombine(VizServer,numP,v);
00651 #ifndef DAGH_NO_MPI
00652 MPI_Barrier(comm_service::comm());
00653 #endif
00654 for (register int l=0; l<=FineLevel(base::GH()); l++) {
00655 Time = CurrentTime(base::GH(),l);
00656 press_cnt = 4;
00657 ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time, l,
00658 press_cnt, base::t[l]);
00659 }
00660 _Interpolation->PointsValuesPar(base::Work(),base::t[0],numP,CloudCtrl[c]._IntPoints,w);
00661 _Interpolation->ArrayCombine(VizServer,numP,w);
00662
00663 for (register int l=0; l<=FineLevel(base::GH()); l++) {
00664 Time = CurrentTime(base::GH(),l);
00665 press_cnt = 1;
00666 ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time, l,
00667 press_cnt, base::t[l]);
00668 }
00669 _Interpolation->PointsValuesPar(base::Work(),base::t[0],numP,CloudCtrl[c]._IntPoints,rho);
00670 _Interpolation->ArrayCombine(VizServer,numP,rho);
00671
00672 fflush(stdout);
00673
00674
00675 if (me == VizServer) {
00676 for (int m=0; m< numP; m++) {
00677 if ( std::sqrt(u[m]*u[m] + v[m]*v[m] + w[m]*w[m] ) > MAXVEL ) {
00678 BBox wb = base::GH().wholebbox();
00679 DCoords wc, min, max, pos;
00680 min = base::GH().worldCoords(wb.lower(),wb.stepsize());
00681 max = base::GH().worldCoords(wb.upper(),wb.stepsize());
00682 wc(0) = CloudCtrl[c]._IntPoints[m](0);
00683 wc(1) = CloudCtrl[c]._IntPoints[m](1);
00684 wc(2) = CloudCtrl[c]._IntPoints[m](2);
00685 if ( inBox(min,max,wc) ) {
00686 fflush(stdout);
00687 ( comm_service::log() << " - - - VELOCITY > MAXVEL at Cloud " << c << " point " << m << std::endl ).flush();
00688 std::cerr << " - - - VELOCITY > MAXVEL at Cloud " << c << " point " << m << std::endl;
00689 assert( std::sqrt(u[m]*u[m] + v[m]*v[m] + w[m]*w[m] ) < MAXVEL );
00690
00691 }
00692 else {
00693 fflush(stdout);
00694 ( comm_service::log() << " - - - Cloud " << c << " point " << m << " EXITS at " << wc << std::endl ).flush();
00695 u[m] = 0.; v[m] = 0.; w[m] = 0.;
00696 }
00697 }
00698 }
00699 std::ofstream outfile;
00700 std::ostream* out;
00701 outfile.open(CloudCtrl[c]._IntName.c_str(), std::ios::out | std::ios::app);
00702 out = new std::ostream(outfile.rdbuf());
00703 *out << base::t[0];
00704 for (int ns=0; ns<numP; ns++) {
00705 *out << " " << CloudCtrl[c]._IntPoints[ns] << " " << rho[ns] << " " << u[ns] << " " << v[ns] << " " << w[ns] << " " << p[ns] << " ";
00706 }
00707 *out << std::endl;
00708 outfile.close();
00709 delete out;
00710
00711 if (nsteps % track_every == 0 ) {
00712 ouputCloud(c,t,rho,u,v,w,p);
00713 }
00714
00715 }
00716
00717 for (int nc=0; nc<numP; nc++) {
00718 CloudCtrl[c]._IntPoints[nc](0) += (u[nc]*dt);
00719 CloudCtrl[c]._IntPoints[nc](1) += (v[nc]*dt);
00720 CloudCtrl[c]._IntPoints[nc](2) += (w[nc]*dt);
00721 CloudCtrl[c].age[nc] += dt;
00722 }
00723
00724 delete [] rho;
00725 delete [] u;
00726 delete [] v;
00727 delete [] w;
00728 delete [] p;
00729 }
00730 }
00731
00732 double modulus(double a, double b) const {
00733 int result = (int) std::floor( a / b );
00734 return a - static_cast<double>( result ) * b;
00735 }
00736
00737 bool inBox (DCoords min, DCoords max, DCoords wc) const {
00738 if ( min(0) <= wc(0) && wc(0) <= max(0)) {
00739 if ( min(1) <= wc(1) && wc(1) <= max(1)) {
00740 if ( min(2) <= wc(2) && wc(2) <= max(2)) {
00741 return 1;
00742 }
00743 }
00744 }
00745 return 0;
00746 }
00747
00748 void ouputCloud(int n, DataType t, DataType* rho, DataType* u, DataType* v, DataType* w, DataType* P) const {
00749 char fname[256];
00750 sprintf(fname,"pointCloud_%02d_%05d.vtk",n,StepsTaken((GridHierarchy&) base::GH(),0) );
00751 std::ofstream ofs;
00752 ofs.open(fname, std::ios::out);
00753 ( comm_service::log() << " *** Writing " << fname << std::endl << std::endl ).flush();
00754
00755 ofs << "# vtk DataFile Version 2.0\n" << "pathLine Test"
00756 << "\nASCII\nDATASET UNSTRUCTURED_GRID\nFIELD FieldData " << 1 << std::endl;
00757 ofs << "TIME 1 1 float\n " << t << "\nPOINTS " << CloudCtrl[n].num_points << " float\n";
00758 for (int i = 0; i < CloudCtrl[n].num_points; i++ ) {
00759 ofs << CloudCtrl[n]._IntPoints[i](0) << " " << CloudCtrl[n]._IntPoints[i](1) << " " << CloudCtrl[n]._IntPoints[i](2) << std::endl;
00760 #ifdef PATH_DEBUG
00761 ( comm_service::log() << t << " pcloud n " << n << " num_points " << CloudCtrl[n].num_points << " " << i << " "
00762 << CloudCtrl[n]._IntPoints[i](0) << " " << CloudCtrl[n]._IntPoints[i](1) << " " << CloudCtrl[n]._IntPoints[i](2) << std::endl ).flush();
00763 #endif
00764 }
00765 ofs << "CELLS " << CloudCtrl[n].num_points << " " << (CloudCtrl[n].num_points)*2 << std::endl;
00766 for (int i = 0; i < (CloudCtrl[n].num_points); i++ ) {
00767 ofs << 1 << " " << i << std::endl;
00768 }
00769 ofs << "CELL_TYPES " << CloudCtrl[n].num_points << std::endl;
00770 for (int i = 0; i < CloudCtrl[n].num_points; i++ )
00771 ofs << 2 << std::endl;
00772 ofs << "POINT_DATA " << CloudCtrl[n].num_points << std::endl;
00773 ofs << "SCALARS rho double\nLOOKUP_TABLE default\n";
00774 for (int i = 0; i < CloudCtrl[n].num_points; i++ ) {
00775 ofs << rho[i] << std::endl;
00776 }
00777 ofs << "SCALARS u double\nLOOKUP_TABLE default\n";
00778 for (int i = 0; i < CloudCtrl[n].num_points; i++ ) {
00779 ofs << u[i] << std::endl;
00780 }
00781 ofs << "SCALARS v double\nLOOKUP_TABLE default\n";
00782 for (int i = 0; i < CloudCtrl[n].num_points; i++ ) {
00783 ofs << v[i] << std::endl;
00784 }
00785 ofs << "SCALARS w double\nLOOKUP_TABLE default\n";
00786 for (int i = 0; i < CloudCtrl[n].num_points; i++ ) {
00787 ofs << w[i] << std::endl;
00788 }
00789 ofs << "SCALARS P double\nLOOKUP_TABLE default\n";
00790 for (int i = 0; i < CloudCtrl[n].num_points; i++ ) {
00791 ofs << P[i] << std::endl;
00792 }
00793 ofs << "SCALARS age double\nLOOKUP_TABLE default\n";
00794 for (int i = 0; i < CloudCtrl[n].num_points; i++ ) {
00795 ofs << CloudCtrl[n].age[i] << std::endl;
00796 }
00797 ofs << "SCALARS cloudNum double\nLOOKUP_TABLE default\n";
00798 for (int i = 0; i < CloudCtrl[n].num_points; i++ ) {
00799 ofs << n << std::endl;
00800 }
00801 ofs.close();
00802 }
00803
00804
00805 virtual bool Restart_(const char* CheckpointFile) {
00806 bool ret = base::Restart_(CheckpointFile);
00807 ( comm_service::log() << " CUSTOM RESTART\n" ).flush();
00808 char fname[256];
00809 sprintf(fname,"%s.cp",CheckpointName.c_str());
00810 std::ifstream ifs(fname, std::ios::in);
00811 sprintf(fname,"viz.cp");
00812 std::ifstream vifs(fname, std::ios::in);
00813 int pos;
00814 vifs.read((char*)&pos,sizeof(int));
00815 vifs.close();
00816 for (register int i=0; i<num_clouds; i++) {
00817 ( comm_service::log() << " pos " << pos << std::endl ).flush();
00818 CloudCtrl[i].Restart(ifs,pos);
00819 pos = ifs.tellg();
00820 }
00821
00822 return ret;
00823 }
00824
00825 virtual void Checkpointing_(const char* CheckpointFile) {
00826 base::Checkpointing_(CheckpointFile);
00827 int me = MY_PROC;
00828 if (me == VizServer) {
00829 ( comm_service::log() << " CUSTOM Viz Checkpoint\n" ).flush();
00830 char fname[256];
00831 sprintf(fname,"%s.cp",CheckpointName.c_str());
00832 std::ofstream ofs(fname, std::ios::app);
00833 int pos = ofs.tellp();
00834 sprintf(fname,"viz.cp");
00835 std::ofstream vofs(fname, std::ios::out);
00836 vofs.write((char*)&pos,sizeof(int));
00837 vofs.close();
00838 for (register int i=0; i<num_clouds; i++) {
00839 CloudCtrl[i].Checkpointing(ofs);
00840 }
00841
00842 }
00843
00844 }
00845
00846 protected:
00847 interpolation_type* _Interpolation;
00848 ControlDevice IntCtrl;
00849
00850 int track_every, num_pathLines, num_clouds;
00851 PointCloudCtrl2 CloudCtrl[maxPointClouds];
00852 };