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 #include "LBMStdProblem.h"
00017 #include "AMRGFMSolver.h"
00018 #include "LBMGFMBoundary.h"
00019 #include "MultiStatCPTLevelSet.h"
00020 #include "AMRGFMInterpolation.h"
00021 #include "Interfaces/SchemeGFMFileOutput.h"
00022 #include "Interfaces/SchemeBrepOutput.h"
00023 #include "LBMD3Q19SupplementalBC.h"
00024
00025
00026 class BoundaryConditionsSpecific : public BoundaryConditionsSupplemental {
00027 typedef BoundaryConditionsSupplemental base;
00028 public:
00029 typedef base::MicroType MicroType;
00030 typedef base::MacroType MacroType;
00031 BoundaryConditionsSpecific(LBMType &lbm) : base(lbm) {}
00032
00033 virtual void SetupData(GridHierarchy* gh, const int& ghosts) {
00034 base::SetupData(gh,ghosts);
00035 DataType l0_p=0.12;
00036 DataType omega = LBM().Omega(LBM().TimeScale());
00037 DataType u_p_avg=(Side(1).Type==7 ? std::sqrt(Side(1).aux[1]*Side(1).aux[1]) : std::sqrt(Side(1).aux[0]*Side(1).aux[0]));
00038 double Re_p=u_p_avg*l0_p/LBM().GasViscosity();
00039 double Re_lbm=(u_p_avg/LBM().VelocityScale())*(l0_p/LBM().LengthScale())/LBM().LatticeViscosity(omega);
00040 ( comm_service::log() << "l0_p=" << l0_p << " u_p_avg=" << u_p_avg
00041 << " Re_p=" << Re_p << " Re_lbm=" << Re_lbm
00042 << "\n omega=" << omega
00043 << " SpeedUp=" << LBM().SpeedUp()
00044 << " u_avg_LBM=" << u_p_avg/LBM().VelocityScale()
00045 << " SpeedRatio=" << u_p_avg/LBM().VelocityScale()*std::sqrt(DataType(3.))
00046 << " LBM().TO()=" << LBM().T0()
00047 << "\n Re_p-Re_lbm=" << Re_p-Re_lbm
00048 << std::endl ).flush();
00049 assert(std::fabs(Re_p-Re_lbm)<DBL_EPSILON*LB_FACTOR*LB_FACTOR*LB_FACTOR);
00050 }
00051
00052 virtual void SetBndry(vec_grid_data_type &fvec, const int& level, const BBox &bb,
00053 const int &dir, const double& time) {
00054
00055 BBox wholebbox = base::GH().wholebbox();
00056 if (dir==LBMType::Right && Side(1).Type==7 && bb.lower(0)==wholebbox.upper(0)+wholebbox.stepsize(0)) {
00057 int Nx = fvec.extents(0), Ny = fvec.extents(1);
00058 int b = fvec.bottom();
00059 assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1)
00060 && fvec.stepsize(2)==bb.stepsize(2));
00061 int mys = std::max(bb.lower(1)/bb.stepsize(1),fvec.lower(1)/fvec.stepsize(1));
00062 int mzs = std::max(bb.lower(2)/bb.stepsize(2),fvec.lower(2)/fvec.stepsize(2));
00063 int mxe = std::min(bb.upper(0)/bb.stepsize(0),fvec.upper(0)/fvec.stepsize(0));
00064 int mye = std::min(bb.upper(1)/bb.stepsize(1),fvec.upper(1)/fvec.stepsize(1));
00065 int mze = std::min(bb.upper(2)/bb.stepsize(2),fvec.upper(2)/fvec.stepsize(2));
00066 MicroType *f = (MicroType *)fvec.databuffer();
00067
00068 DataType rho=Side(1).aux[0];
00069 DataType a0=Side(1).aux[1];
00070 DataType a1=Side(1).aux[2];
00071 DataType a2=Side(1).aux[3];
00072
00073 if (Side(1).Scaling==LBMType::Physical) {
00074 rho /= LBM().DensityScale();
00075 a0 /= LBM().VelocityScale();
00076 a1 /= LBM().VelocityScale();
00077 a2 /= LBM().VelocityScale();
00078 }
00079
00080 MacroType q;
00081 q(0) = rho;
00082 q(1) = a0;
00083 q(2) = a1;
00084 q(3) = a2;
00085
00086 for (register int k=mzs; k<=mze; k++)
00087 for (register int j=mys; j<=mye; j++)
00088 f[b+LBM().idx(mxe,j,k,Nx,Ny)] = LBM().Equilibrium(q);
00089 }
00090 else
00091 base::SetBndry(fvec,level,bb,dir,time);
00092
00093 }
00094 };
00095
00096
00097 class SolverSpecific :
00098 public AMRGFMSolver<MicroType,FixupType,FlagType,DIM> {
00099 typedef AMRGFMSolver<MicroType,FixupType,FlagType,DIM> base;
00100 typedef SchemeBrepOutput<LBMType,DIM> brep_output_type;
00101 typedef AMRGFMInterpolation<MicroType,FixupType,FlagType,DIM> interpolation_type;
00102 typedef MicroType::InternalDataType DataType;
00103 typedef brep_output_type::geom_point_type geom_point_type;
00104 typedef SchemeGFMFileOutput<LBMType,FixupType,FlagType,DIM> output_type;
00105 typedef MultiStatCPTLevelSet<DataType,DIM> cpt_type;
00106 public:
00107 SolverSpecific(IntegratorSpecific& integ,
00108 base::initial_condition_type& init,
00109 base::boundary_conditions_type& bc) : base(integ, init, bc) {
00110 SetLevelTransfer(new LBMF77LevelTransfer<LBMType,DIM>(integ.Scheme(), f_prolong, f_restrict));
00111 SetFileOutput(new output_type(*this,integ.Scheme()));
00112 SetFixup(new FixupSpecific(integ.Scheme()));
00113 SetFlagging(new FlaggingSpecific(*this,integ.Scheme()));
00114 AddGFM(new GhostFluidMethod<MicroType,DIM>(new LBMGFMBoundary<LBMType,DIM>(integ.LBM()),
00115 new cpt_type()));
00116 _BrepOutput = new brep_output_type(new interpolation_type(*this),integ.Scheme());
00117 IntegrateEvery = 1;
00118 }
00119
00120 ~SolverSpecific() {
00121 DeleteGFM(_GFM[0]);
00122 delete _LevelTransfer;
00123 delete _Flagging;
00124 delete _Fixup;
00125 delete _FileOutput;
00126 delete _BrepOutput;
00127 }
00128
00129 virtual void SetupData() {
00130 base::SetupData();
00131 base::NAMRTimeSteps = 1;
00132 base::AdaptBndTimeInterpolate = 0;
00133 base::Step[0].LastTime = 1.e37;
00134 base::Step[0].VariableTimeStepping = -1;
00135 base::Step[0].dtv[0] = ((LBMIntegrator<LBMType,DIM> &)Integrator_()).LBM().TimeScale();
00136 if (_BrepOutput) _BrepOutput->SetupData(base::PGH(), base::NGhosts());
00137 }
00138
00139 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00140 base::register_at(Ctrl,prefix);
00141 if (_BrepOutput) _BrepOutput->register_at(base::LocCtrl,prefix);
00142 RegisterAt(base::LocCtrl,"OutputDragEvery",IntegrateEvery);
00143 }
00144
00145 virtual void update() {
00146 base::update();
00147 if (_BrepOutput) _BrepOutput->update();
00148 }
00149
00150 virtual void Advance(double& t, double& dt) {
00151 base::Advance(t,dt);
00152 if (IntegrateEvery>0)
00153 if (StepsTaken(base::GH(),0)%IntegrateEvery == 0) {
00154 int num_vertices=0, num_connections=0;
00155 const DataType* vertices=0;
00156 const int* connections=0;
00157 if (!((cpt_type *)base::GFM(0).GetLevelSetP())->GetNthBrep(0,num_vertices,vertices,num_connections,connections))
00158 return;
00159
00160
00161 for (register int l=0; l<=FineLevel(base::GH()); l++) {
00162 int Time = CurrentTime(base::GH(),l);
00163 int press_cnt = base::Dim()+4;
00164 ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time, l,
00165 press_cnt, base::t[l]);
00166 }
00167
00168 geom_point_type moment_origin = ((cpt_type *)base::GFM(0).GetLevelSetP())->morigin[0];
00169 geom_point_type force = _BrepOutput->IntegrateNormal(base::Work(),num_vertices,vertices,num_connections,connections);
00170 geom_point_type moment = _BrepOutput->IntegrateMoment(base::Work(),num_vertices,vertices,num_connections,connections,moment_origin);
00171
00172 if (MY_PROC==VizServer) {
00173 std::ofstream outfile;
00174 std::ostream* out;
00175 std::string fname = "drag.dat";
00176 outfile.open(fname.c_str(), std::ios::out | std::ios::app);
00177 out = new std::ostream(outfile.rdbuf());
00178 *out << StepsTaken(base::GH(),0) << " " << t << " " << force << " " << moment << std::endl;
00179 outfile.close();
00180 delete out;
00181 }
00182 }
00183 }
00184
00185 virtual void Output() {
00186 if (CurrentTime(base::GH(),0) == base::LastOutputTime()) return;
00187
00188 int num_vertices=0, num_edges=0;
00189 const DataType* vertices=0;
00190 const int* edges=0;
00191
00192 if (!((MultiStatCPTLevelSet<DataType,DIM> *)base::GFM(0).GetLevelSetP())->GetBrep(num_vertices,vertices,num_edges,edges))
00193 return;
00194
00195 _BrepOutput->WriteOut(base::U(),num_vertices,vertices,num_edges,edges);
00196
00197 base::Output();
00198 }
00199
00200 protected:
00201 brep_output_type* _BrepOutput;
00202 int IntegrateEvery;
00203 };
00204
00205 #endif