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 BBox wholebbox = base::GH().wholebbox();
00055 if (dir==LBMType::Front && Side(4).Type==20 && bb.upper(2)+bb.stepsize(2)==wholebbox.lower(2) &&
00056 Side(4).NAux==4) {
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 mxs = std::max(bb.lower(0)/bb.stepsize(0),fvec.lower(0)/fvec.stepsize(0));
00062 int mys = std::max(bb.lower(1)/bb.stepsize(1),fvec.lower(1)/fvec.stepsize(1));
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 DCoords lbcorner = base::GH().worldCoords(fvec.lower(), fvec.stepsize());
00069 DCoords dx = base::GH().worldStep(fvec.stepsize());
00070 DataType xs=Side(4).aux[0]; DataType xe=Side(4).aux[1];
00071 DataType ys=Side(4).aux[2]; DataType ye=Side(4).aux[3];
00072
00073 for (register int j=mys; j<=mye; j++)
00074 for (register int i=mxs; i<=mxe; i++) {
00075 double xc = (i+0.5)*dx(0)+lbcorner(0);
00076 double yc = (j+0.5)*dx(1)+lbcorner(1);
00077
00078 if (xc>=xs && xc<=xe && yc>=ys && yc<=ye) {
00079
00080 if (i>mxs) f[b+base::LBM().idx(i,j,mze,Nx,Ny)](18) = f[b+base::LBM().idx(i-1,j,mze+1,Nx,Ny)](17);
00081 f[b+base::LBM().idx(i,j,mze,Nx,Ny)](5) = f[b+base::LBM().idx(i,j,mze+1,Nx,Ny)](6);
00082 if (i<mxe) f[b+base::LBM().idx(i,j,mze,Nx,Ny)](15) = f[b+base::LBM().idx(i+1,j,mze+1,Nx,Ny)](16);
00083
00084 if (j>mys) f[b+base::LBM().idx(i,j,mze,Nx,Ny)](14) = f[b+base::LBM().idx(i,j-1,mze+1,Nx,Ny)](13);
00085 if (j<mye) f[b+base::LBM().idx(i,j,mze,Nx,Ny)](11) = f[b+base::LBM().idx(i,j+1,mze+1,Nx,Ny)](12);
00086 }
00087 else {
00088
00089 if (i>mxs) f[b+base::LBM().idx(i,j,mze,Nx,Ny)](15) = f[b+base::LBM().idx(i-1,j,mze+1,Nx,Ny)](17);
00090 f[b+base::LBM().idx(i,j,mze,Nx,Ny)](5) = f[b+base::LBM().idx(i,j,mze+1,Nx,Ny)](6);
00091 if (i<mxe) f[b+base::LBM().idx(i,j,mze,Nx,Ny)](18) = f[b+base::LBM().idx(i+1,j,mze+1,Nx,Ny)](16);
00092
00093 if (j>mys) f[b+base::LBM().idx(i,j,mze,Nx,Ny)](11) = f[b+base::LBM().idx(i,j-1,mze+1,Nx,Ny)](13);
00094 if (j<mye) f[b+base::LBM().idx(i,j,mze,Nx,Ny)](14) = f[b+base::LBM().idx(i,j+1,mze+1,Nx,Ny)](12);
00095 }
00096 }
00097 }
00098 else
00099 base::SetBndry(fvec,level,bb,dir,time);
00100 }
00101 };
00102
00103
00104 class SolverSpecific :
00105 public AMRGFMSolver<MicroType,FixupType,FlagType,DIM> {
00106 typedef AMRGFMSolver<MicroType,FixupType,FlagType,DIM> base;
00107 typedef SchemeBrepOutput<LBMType,DIM> brep_output_type;
00108 typedef AMRGFMInterpolation<MicroType,FixupType,FlagType,DIM> interpolation_type;
00109 typedef MicroType::InternalDataType DataType;
00110 typedef brep_output_type::geom_point_type geom_point_type;
00111 typedef SchemeGFMFileOutput<LBMType,FixupType,FlagType,DIM> output_type;
00112 typedef MultiStatCPTLevelSet<DataType,DIM> cpt_type;
00113 public:
00114 SolverSpecific(IntegratorSpecific& integ,
00115 base::initial_condition_type& init,
00116 base::boundary_conditions_type& bc) : base(integ, init, bc) {
00117 SetLevelTransfer(new LBMF77LevelTransfer<LBMType,DIM>(integ.Scheme(), f_prolong, f_restrict));
00118 SetFileOutput(new output_type(*this,integ.Scheme()));
00119 SetFixup(new FixupSpecific(integ.Scheme()));
00120 SetFlagging(new FlaggingSpecific(*this,integ.Scheme()));
00121 AddGFM(new GhostFluidMethod<MicroType,DIM>(new LBMGFMBoundary<LBMType,DIM>(integ.LBM()),
00122 new cpt_type()));
00123 _BrepOutput = new brep_output_type(new interpolation_type(*this),integ.Scheme());
00124 IntegrateEvery = 0;
00125 moment_origin = 0.;
00126 qref = aref = href = 1.;
00127 }
00128
00129 ~SolverSpecific() {
00130 DeleteGFM(_GFM[0]);
00131 delete _LevelTransfer;
00132 delete _Flagging;
00133 delete _Fixup;
00134 delete _FileOutput;
00135 delete _BrepOutput;
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 if (_BrepOutput) _BrepOutput->SetupData(base::PGH(), base::NGhosts());
00146
00147
00148 DataType *aux = ((BoundaryConditionsSpecific &)BoundaryConditions_()).Side(0).aux;
00149 if (((BoundaryConditionsSpecific &)BoundaryConditions_()).Side(0).Type==7) qref=0.5*aux[0]*(aux[1]*aux[1]);
00150 if (((BoundaryConditionsSpecific &)BoundaryConditions_()).Side(0).Type==3)
00151 qref=0.5*((LBMIntegrator<LBMType,DIM> &)Integrator_()).LBM().GasDensity()*(aux[0]*aux[0]);
00152 moment_origin = ((cpt_type *)base::GFM(0).GetLevelSetP())->ForwardTransformN(0,moment_origin);
00153 ( comm_service::log() << "qref=" << qref << " moment_origin=" << moment_origin << std::endl ).flush();
00154 }
00155
00156 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00157 base::register_at(Ctrl,prefix);
00158 if (_BrepOutput) _BrepOutput->register_at(base::LocCtrl,prefix);
00159 RegisterAt(base::LocCtrl,"OutputDragEvery",IntegrateEvery);
00160 RegisterAt(base::LocCtrl,"Aref",aref);
00161 RegisterAt(base::LocCtrl,"Href",href);
00162 char VariableName[32];
00163 for (int d=0; d<base::Dim(); d++) {
00164 std::sprintf(VariableName,"MOrigin(%d)",d+1);
00165 RegisterAt(base::LocCtrl,VariableName,moment_origin[d]);
00166 }
00167 }
00168
00169 virtual void update() {
00170 base::update();
00171 if (_BrepOutput) _BrepOutput->update();
00172 }
00173
00174 virtual void Advance(double& t, double& dt) {
00175 base::Advance(t,dt);
00176 if (IntegrateEvery>0)
00177 if (StepsTaken(base::GH(),0)%IntegrateEvery == 0) {
00178 int num_vertices=0, num_connections=0;
00179 const DataType* vertices=0;
00180 const int* connections=0;
00181 if (!((cpt_type *)base::GFM(0).GetLevelSetP())->GetNthBrep(0,num_vertices,vertices,num_connections,connections))
00182 return;
00183
00184
00185 for (register int l=0; l<=FineLevel(base::GH()); l++) {
00186 int Time = CurrentTime(base::GH(),l);
00187 int press_cnt = base::Dim()+4;
00188 ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time, l,
00189 press_cnt, base::t[l]);
00190 }
00191
00192 geom_point_type force = _BrepOutput->IntegrateNormal(base::Work(),num_vertices,vertices,num_connections,connections);
00193 geom_point_type moment = _BrepOutput->IntegrateMoment(base::Work(),num_vertices,vertices,num_connections,connections,moment_origin);
00194
00195 if (MY_PROC==VizServer) {
00196 std::ofstream outfile;
00197 std::ostream* out;
00198 std::string fname = "drag.dat";
00199 outfile.open(fname.c_str(), std::ios::out | std::ios::app);
00200 out = new std::ostream(outfile.rdbuf());
00201
00202
00203 *out << StepsTaken(base::GH(),0) << " " << t << " " << force << " " << moment << " ";
00204
00205
00206 DataType scale = ((cpt_type *)base::GFM(0).GetLevelSetP())->scale[0];
00207 force = ((cpt_type *)base::GFM(0).GetLevelSetP())->InverseTransformN(0,force)*scale;
00208 moment = ((cpt_type *)base::GFM(0).GetLevelSetP())->InverseTransformN(0,moment)*scale;
00209 *out << force/(aref*qref) << " " << moment/(aref*qref*href) << std::endl;
00210
00211 outfile.close();
00212 delete out;
00213 }
00214 }
00215 }
00216
00217 virtual void Output() {
00218 if (CurrentTime(base::GH(),0) == base::LastOutputTime()) return;
00219
00220 int num_vertices=0, num_edges=0;
00221 const DataType* vertices=0;
00222 const int* edges=0;
00223
00224 if (!((MultiStatCPTLevelSet<DataType,DIM> *)base::GFM(0).GetLevelSetP())->GetBrep(num_vertices,vertices,num_edges,edges))
00225 return;
00226
00227 _BrepOutput->WriteOut(base::U(),num_vertices,vertices,num_edges,edges);
00228
00229 base::Output();
00230 }
00231
00232 protected:
00233 brep_output_type* _BrepOutput;
00234 int IntegrateEvery;
00235 geom_point_type moment_origin;
00236 DataType qref, aref, href;
00237 };
00238
00239 #endif