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_LBMSCHEME
00015 #define OWN_AMRSOLVER
00016 #include "LBMStdProblem.h"
00017 #include "AMRGFMSolver.h"
00018 #include "LBMGFMBoundary.h"
00019 #include "StatCPTLevelSet.h"
00020 #include "AMRGFMInterpolation.h"
00021 #include "Interfaces/SchemeGFMFileOutput.h"
00022 #include "Interfaces/SchemeBrepOutput.h"
00023
00024 DataType xs[DIM], 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 << "D3Q19: 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 GFMBoundarySpecific : public LBMGFMBoundary<LBMType,DIM> {
00058 typedef LBMGFMBoundary<LBMType,DIM> base;
00059 public:
00060 GFMBoundarySpecific(LBMType &lbm) : base(lbm) {}
00061
00062 virtual void SetBndry(vec_grid_data_type& gdu, const int& nc, const int* idx,
00063 const MicroType* u, const point_type* xc, const DataType* distance,
00064 const point_type* normal, DataType* aux, const int naux, const double t)
00065 { LBM().GFMBCStandard(gdu,base::Type,nc,idx,u,xc,distance,normal,aux,naux,LBMType::Lattice); }
00066
00067 virtual void SetBndryAux(vec_grid_data_type& gdu, grid_data_type& gdphi,
00068 const MicroType* u, DataType* aux, const int& Level,
00069 double t, const int& nc, const int* idx,
00070 const point_type* xc, const DataType* distance,
00071 const point_type* normal) {
00072 for (register int i=0; i<nc; i++) {
00073 DataType x = xc[i](0)-distance[i]*normal[i](0)-xs[0];
00074 DataType y = xc[i](1)-distance[i]*normal[i](1)-xs[1];
00075 DataType rs = std::sqrt(x*x+y*y);
00076 DataType alpha = std::atan2(y,x);
00077 aux[i*base::NAux()] = 0.;
00078 aux[i*base::NAux()+1] = 0.;
00079 aux[i*base::NAux()+2] = 0.;
00080 if (Tr!=0.) {
00081 DataType vr = 2.*pi/Tr*rs;
00082 aux[i*base::NAux()] -= vr*std::sin(alpha);
00083 aux[i*base::NAux()+1] += vr*std::cos(alpha);
00084 }
00085 }
00086 }
00087 };
00088
00089
00090 class GFMLevelSetSpecific : public StatCPTLevelSet<DataType,DIM> {
00091 typedef StatCPTLevelSet<DataType,DIM> base;
00092
00093 public:
00094 typedef base::grid_fct_type grid_fct_type;
00095 typedef base::grid_data_type grid_data_type;
00096
00097 GFMLevelSetSpecific() : scale(1.), num_vertices(0), initial_vertices(0) {}
00098 virtual ~GFMLevelSetSpecific() {
00099 if (initial_vertices) delete [] initial_vertices;
00100 }
00101
00102 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00103 base::register_at(Ctrl,prefix);
00104 RegisterAt(base::LocCtrl,"Scale",scale);
00105 }
00106
00107 virtual void SetupData(GridHierarchy* gh, const int& ghosts) {
00108 base::SetupData(gh, ghosts);
00109 base::_Stationary = 0;
00110 int num_connections;
00111 const int* connections;
00112 const DataType* ver;
00113 base::GetBrep(num_vertices, ver, num_connections, connections);
00114 point_type* vertices = (point_type*)(ver);
00115 initial_vertices = new point_type[ num_vertices ];
00116 for (int i=0; i<num_vertices; i++) {
00117 vertices[i] *= scale;
00118 initial_vertices[i] = vertices[i];
00119 }
00120 base::SetBrep(num_vertices, ver, num_connections, connections);
00121 }
00122
00123 virtual void SetPhi(grid_fct_type& phi, const int Time, const int Level, double t) {
00124 if (Tr!=0.) {
00125 int num_connections;
00126 const int* connections;
00127 const DataType* ver;
00128 base::GetBrep(num_vertices, ver, num_connections, connections);
00129 point_type* vertices = (point_type*)(ver);
00130
00131 for (int i=0; i<num_vertices; i++) {
00132 DataType x = initial_vertices[i](0)-xs[0];
00133 DataType y = initial_vertices[i](1)-xs[1];
00134 DataType rs = std::sqrt(x*x+y*y);
00135 DataType alpha = std::atan2(y,x);
00136
00137 vertices[i](0) = xs[0]+rs*std::cos(2.*pi/Tr*t+alpha);
00138 vertices[i](1) = xs[1]+rs*std::sin(2.*pi/Tr*t+alpha);
00139 vertices[i](2) = initial_vertices[i](2);
00140 }
00141 base::SetBrep(num_vertices, ver, num_connections, connections);
00142 }
00143 base::SetPhi(phi, Time, Level, t);
00144 }
00145
00146 protected:
00147 DataType scale;
00148 int num_vertices;
00149 point_type* initial_vertices;
00150 };
00151
00152
00153 class GhostFluidMethodSpecific : public GhostFluidMethod<MicroType,DIM> {
00154 typedef GhostFluidMethod<MicroType,DIM> base;
00155 public:
00156 GhostFluidMethodSpecific(boundary_type* boundary, levelset_type* levelset) :
00157 base(boundary,levelset) {
00158 xs[0]=0.; xs[1]=0.; xs[2]=0.; Tr=0.;
00159 }
00160
00161 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00162 base::register_at(Ctrl,prefix);
00163 RegisterAt(base::LocCtrl,"Center(1)",xs[0]);
00164 RegisterAt(base::LocCtrl,"Center(2)",xs[1]);
00165 RegisterAt(base::LocCtrl,"Center(3)",xs[2]);
00166 RegisterAt(base::LocCtrl,"RotationDuration",Tr);
00167 }
00168 };
00169
00170
00171 class SolverSpecific :
00172 public AMRGFMSolver<MicroType,FixupType,FlagType,DIM> {
00173 typedef AMRGFMSolver<MicroType,FixupType,FlagType,DIM> base;
00174 typedef SchemeBrepOutput<LBMType,DIM> brep_output_type;
00175 typedef AMRGFMInterpolation<MicroType,FixupType,FlagType,DIM> interpolation_type;
00176 typedef MicroType::InternalDataType DataType;
00177 public:
00178 SolverSpecific(IntegratorSpecific& integ,
00179 base::initial_condition_type& init,
00180 base::boundary_conditions_type& bc) : base(integ, init, bc) {
00181 SetLevelTransfer(new LBMF77LevelTransfer<LBMType,DIM>(integ.Scheme(), f_prolong, f_restrict));
00182 SetFileOutput(new SchemeGFMFileOutput<LBMType,FixupType,FlagType,DIM>(*this,integ.Scheme()));
00183 SetFixup(new FixupSpecific(integ.Scheme()));
00184 SetFlagging(new FlaggingSpecific(*this,integ.Scheme()));
00185 AddGFM(new GhostFluidMethodSpecific(new GFMBoundarySpecific(integ.LBM()),
00186 new GFMLevelSetSpecific()));
00187 _BrepOutput = new brep_output_type(new interpolation_type(*this),integ.Scheme());
00188 }
00189
00190 ~SolverSpecific() {
00191 DeleteGFM(_GFM[0]);
00192 delete _LevelTransfer;
00193 delete _Flagging;
00194 delete _Fixup;
00195 delete _FileOutput;
00196 delete _BrepOutput;
00197 }
00198
00199 virtual void SetupData() {
00200 base::SetupData();
00201 base::NAMRTimeSteps = 1;
00202 base::AdaptBndTimeInterpolate = 0;
00203 base::Step[0].LastTime = 1.e37;
00204 base::Step[0].VariableTimeStepping = -1;
00205 base::Step[0].dtv[0] = ((LBMIntegrator<LBMType,DIM> &)Integrator_()).LBM().TimeScale();
00206 if (_BrepOutput) _BrepOutput->SetupData(base::PGH(), base::NGhosts());
00207 }
00208
00209 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00210 base::register_at(Ctrl,prefix);
00211 if (_BrepOutput) _BrepOutput->register_at(base::LocCtrl,prefix);
00212 }
00213
00214 virtual void update() {
00215 base::update();
00216 if (_BrepOutput) _BrepOutput->update();
00217 }
00218
00219 virtual void Output() {
00220 if (CurrentTime(base::GH(),0) == base::LastOutputTime()) return;
00221
00222 int num_vertices=0, num_edges=0;
00223 const DataType* vertices=0;
00224 const int* edges=0;
00225
00226 if (!((CPTLevelSet<DataType,DIM> *)base::GFM(0).GetLevelSetP())->GetBrep(num_vertices,vertices,num_edges,edges))
00227 return;
00228
00229 _BrepOutput->WriteOut(base::U(),num_vertices,vertices,num_edges,edges);
00230
00231 base::Output();
00232 }
00233
00234 protected:
00235 brep_output_type* _BrepOutput;
00236 };
00237
00238 #endif