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_LBMSCHEME
00015 #define OWN_INITIALCONDITION
00016 #define OWN_AMRSOLVER
00017 #include "LBMStdProblem.h"
00018 #include "AMRGFMSolver.h"
00019 #include "LBMGFMBoundary.h"
00020 #include "Interfaces/SchemeGFMFileOutput.h"
00021 #include "StatCPTLevelSet.h"
00022
00023 DataType rho0=1.;
00024 DataType xs[2], Tr;
00025 DataType pi=4.0*std::atan(1.);
00026
00027 class LBMSpecific : public LBMType {
00028 typedef LBMType base;
00029 public:
00030 LBMSpecific() : base(), omega(1.) {}
00031
00032 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00033 base::register_at(Ctrl,prefix);
00034 RegisterAt(base::LocCtrl,"Omega",omega);
00035 }
00036
00037 virtual void SetupData(GridHierarchy* gh, const int& ghosts) {
00038 base::SetupData(gh,ghosts);
00039 if (base::LengthScale()!=DataType(1.) || base::TimeScale()!=DataType(1.)) {
00040 std::cerr << "dx(0) and dt(0) need to be equal to 1.0!" << std::endl;
00041 std::exit(-1);
00042 }
00043 if (omega<0.5 || omega>2.) {
00044 std::cerr << "0.5 <= omega <= 2 required." << std::endl;
00045 std::exit(-1);
00046 }
00047 base::SetGas(1.,base::LatticeViscosity(omega),base::LatticeSpeedOfSound());
00048 std::cout << "D2Q9: Gas_rho=" << base::GasDensity() << " Gas_Cs=" << base::GasSpeedofSound()
00049 << " Gas_nu=" << base::GasViscosity() << std::endl;
00050 }
00051
00052 private:
00053 DataType omega;
00054 };
00055
00056
00057 class InitialConditionSpecific : public LBMInitialCondition<LBMType,DIM> {
00058 typedef LBMInitialCondition<LBMType,DIM> base;
00059 public:
00060 InitialConditionSpecific(LBMType &lbm) : base(lbm) {}
00061
00062 virtual void SetGrid(vec_grid_data_type& fvec, grid_data_type& workvec, const int& level) {
00063 DataType q[3] = { rho0, 0., 0. };
00064 LBM().ICStandard(fvec,LBMType::ConstantMacro,q,3,LBMType::Lattice);
00065 }
00066 };
00067
00068
00069
00070 class GFMBoundarySpecific : public LBMGFMBoundary<LBMType,DIM> {
00071 typedef LBMGFMBoundary<LBMType,DIM> base;
00072 public:
00073 GFMBoundarySpecific(LBMType &lbm) : base(lbm) {}
00074
00075 virtual void SetBndry(vec_grid_data_type& gdu, const int& nc, const int* idx,
00076 const MicroType* u, const point_type* xc, const DataType* distance,
00077 const point_type* normal, DataType* aux, const int naux, const double t)
00078 { LBM().GFMBCStandard(gdu,base::Type,nc,idx,u,xc,distance,normal,aux,naux,LBMType::Lattice); }
00079
00080 virtual void SetBndryAux(vec_grid_data_type& gdu, grid_data_type& gdphi,
00081 const MicroType* u, DataType* aux, const int& Level,
00082 double t, const int& nc, const int* idx,
00083 const point_type* xc, const DataType* distance,
00084 const point_type* normal) {
00085 for (register int i=0; i<nc; i++) {
00086 DataType x = xc[i](0)-distance[i]*normal[i](0)-xs[0];
00087 DataType y = xc[i](1)-distance[i]*normal[i](1)-xs[1];
00088 DataType rs = std::sqrt(x*x+y*y);
00089 DataType alpha = std::atan2(y,x);
00090 aux[i*base::NAux()] = 0.;
00091 aux[i*base::NAux()+1] = 0.;
00092 if (Tr!=0.) {
00093 DataType vr = 2.*pi/Tr*rs;
00094 aux[i*base::NAux()] -= vr*std::sin(alpha);
00095 aux[i*base::NAux()+1] += vr*std::cos(alpha);
00096 }
00097 }
00098 }
00099 };
00100
00101
00102 class GFMLevelSetSpecific : public StatCPTLevelSet<DataType,DIM> {
00103 typedef StatCPTLevelSet<DataType,DIM> base;
00104
00105 public:
00106 typedef base::grid_fct_type grid_fct_type;
00107 typedef base::grid_data_type grid_data_type;
00108
00109 GFMLevelSetSpecific() : scale(1.), num_vertices(0), initial_vertices(0) {}
00110 virtual ~GFMLevelSetSpecific() {
00111 if (initial_vertices) delete [] initial_vertices;
00112 }
00113
00114 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00115 base::register_at(Ctrl,prefix);
00116 RegisterAt(base::LocCtrl,"Scale",scale);
00117 }
00118
00119 virtual void SetupData(GridHierarchy* gh, const int& ghosts) {
00120 base::SetupData(gh, ghosts);
00121 base::_Stationary = 0;
00122 int num_connections;
00123 const int* connections;
00124 const DataType* ver;
00125 base::GetBrep(num_vertices, ver, num_connections, connections);
00126 point_type* vertices = (point_type*)(ver);
00127 initial_vertices = new point_type[ num_vertices ];
00128 for (int i=0; i<num_vertices; i++) {
00129 vertices[i] *= scale;
00130 initial_vertices[i] = vertices[i];
00131 }
00132 base::SetBrep(num_vertices, ver, num_connections, connections);
00133 }
00134
00135 virtual void SetPhi(grid_fct_type& phi, const int Time, const int Level, double t) {
00136 if (Tr!=0.) {
00137 int num_connections;
00138 const int* connections;
00139 const DataType* ver;
00140 base::GetBrep(num_vertices, ver, num_connections, connections);
00141 point_type* vertices = (point_type*)(ver);
00142
00143 for (int i=0; i<num_vertices; i++) {
00144 DataType x = initial_vertices[i](0)-xs[0];
00145 DataType y = initial_vertices[i](1)-xs[1];
00146 DataType rs = std::sqrt(x*x+y*y);
00147 DataType alpha = std::atan2(y,x);
00148
00149 vertices[i](0) = xs[0]+rs*std::cos(2.*pi/Tr*t+alpha);
00150 vertices[i](1) = xs[1]+rs*std::sin(2.*pi/Tr*t+alpha);
00151 }
00152 base::SetBrep(num_vertices, ver, num_connections, connections);
00153 }
00154 base::SetPhi(phi, Time, Level, t);
00155 }
00156
00157 protected:
00158 DataType scale;
00159 int num_vertices;
00160 point_type* initial_vertices;
00161 };
00162
00163
00164 class GhostFluidMethodSpecific : public GhostFluidMethod<MicroType,DIM> {
00165 typedef GhostFluidMethod<MicroType,DIM> base;
00166 public:
00167 GhostFluidMethodSpecific(boundary_type* boundary, levelset_type* levelset) :
00168 base(boundary,levelset) {
00169 xs[0]=0.; xs[1]=0.; Tr=0.;
00170 }
00171
00172 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00173 base::register_at(Ctrl,prefix);
00174 RegisterAt(base::LocCtrl,"Center(1)",xs[0]);
00175 RegisterAt(base::LocCtrl,"Center(2)",xs[1]);
00176 RegisterAt(base::LocCtrl,"RotationDuration",Tr);
00177 }
00178 };
00179
00180
00181 class SolverSpecific :
00182 public AMRGFMSolver<MicroType,FixupType,FlagType,DIM> {
00183 typedef AMRGFMSolver<MicroType,FixupType,FlagType,DIM> base;
00184 public:
00185 SolverSpecific(IntegratorSpecific& integ,
00186 base::initial_condition_type& init,
00187 base::boundary_conditions_type& bc) : base(integ, init, bc) {
00188 SetLevelTransfer(new LBMF77LevelTransfer<LBMType,DIM>(integ.Scheme(), f_prolong, f_restrict));
00189 SetFileOutput(new SchemeGFMFileOutput<LBMType,FixupType,FlagType,DIM>(*this,integ.Scheme()));
00190 SetFixup(new FixupSpecific(integ.Scheme()));
00191 SetFlagging(new FlaggingSpecific(*this,integ.Scheme()));
00192 AddGFM(new GhostFluidMethodSpecific(new GFMBoundarySpecific(integ.LBM()),
00193 new GFMLevelSetSpecific()));
00194 }
00195
00196 ~SolverSpecific() {
00197 DeleteGFM(_GFM[0]);
00198 delete _LevelTransfer;
00199 delete _Flagging;
00200 delete _Fixup;
00201 delete _FileOutput;
00202 }
00203
00204 virtual void SetupData() {
00205 base::SetupData();
00206 base::NAMRTimeSteps = 1;
00207 base::AdaptBndTimeInterpolate = 0;
00208 base::Step[0].LastTime = 1.e37;
00209 base::Step[0].VariableTimeStepping = -1;
00210 base::Step[0].dtv[0] = ((LBMIntegrator<LBMType,DIM> &)Integrator_()).LBM().TimeScale();
00211 }
00212 };
00213
00214 #endif