00001
00002
00003
00004
00005
00006 #ifndef AMROC_LBM_FIXUP_H
00007 #define AMROC_LBM_FIXUP_H
00008
00016 #include "AMRFixup.h"
00017 #include "LBMFixupOps.h"
00018
00019 template <class LBMType, class FixupType, int dim>
00020 class LBMFixup : public AMRFixup<typename LBMType::MicroType,FixupType,dim>,
00021 public LBMFixupOps<LBMType,dim> {
00022 typedef AMRFixup<typename LBMType::MicroType,FixupType,dim> base;
00023 typedef AMRBase<typename LBMType::MicroType,dim> amrbase;
00024 typedef LBMFixupOps<LBMType,dim> lbmfixupops;
00025
00026 public:
00027 typedef typename LBMType::MicroType MicroType;
00028 typedef typename LBMType::MacroType MacroType;
00029 typedef typename LBMType::DataType DataType;
00030
00031 typedef typename base::vec_grid_data_type vec_grid_data_type;
00032
00033 LBMFixup(LBMType &lbm) : _lbm(lbm), _RescaleNonEq(0) {
00034 for (int d=0; d<2*dim; d++)
00035 indices[d] = (int *) 0;
00036 }
00037 virtual ~LBMFixup() {
00038 for (int d=0; d<2*dim; d++)
00039 if (indices[d]) delete [] indices[d];
00040 }
00041
00042 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00043 base::LocCtrl = Ctrl.getSubDevice(prefix+"Fixup");
00044 RegisterAt(base::LocCtrl,"RescaleNonEq",_RescaleNonEq);
00045 }
00046 virtual void register_at(ControlDevice& Ctrl) { register_at(Ctrl,""); }
00047
00048 virtual void SetupData(GridHierarchy* gh, const int& ghosts) {
00049 amrbase::SetupData(gh, ghosts);
00050 for (int d=0; d<2*dim; d++) {
00051 indices[d] = new int[LBM().NMicroVar()];
00052 Nindices[d] = LBM().OutgoingIndices(d,indices[d]);
00053 for (register int n=Nindices[d]; n<LBM().NMicroVar(); n++)
00054 indices[d][n]=-1;
00055 }
00056 }
00057
00058 virtual void SaveFluxes(const int Time, const int Level, const int c,
00059 vec_grid_data_type* flux[],
00060 const double dt, const int& mpass) {}
00061
00062 virtual void AddIntercellFluxes(const int Time, const int Level, const int c,
00063 vec_grid_data_type* flux[], const double dt,
00064 const int& mpass) {}
00065
00066 virtual void Correction(const int Time, const int WTime, const int Level,
00067 const double t, const double dt) {
00068
00069 int TimeF = CurrentTime(base::GH(),Level+1);
00070 int TimeC = CurrentTime(base::GH(),Level);
00071 const int refinefactor = RefineFactor(base::GH(),Level);
00072 DataType fac = refinefactor*LBM().Omega(dt/refinefactor)/LBM().Omega(dt);
00073
00074 forall (base::U(),TimeF,Level+1,c)
00075 for (int d=0; d<2*dim; d++) {
00076 if (base::ValidSide(TimeF,Level+1,c,d)) {
00077 BBox bb = base::U()(TimeF,Level+1,c).bbox();
00078 BBox from(CorrectionBBox(coarsen(base::U().interiorbbox(TimeF,Level+1,c),refinefactor),d));
00079 if (d%2 == 0) from.growupper(d/2,1);
00080 else from.growlower(d/2,-1);
00081 BBox to(base::FixupBBox(coarsen(base::U().interiorbbox(TimeF,Level+1,c),refinefactor),d));
00082 vec_grid_data_type uh(from);
00083 lbmfixupops::average_indices(uh,base::U()(TimeF,Level+1,c),from,indices[d],Nindices[d]);
00084 LBM().ReverseStream(uh,to,d);
00085 for (register int j=0; j<base::U().parents(TimeF,Level+1,c); j++) {
00086 BBox wto(coarsen(base::U().parentbox(TimeF,Level+1,c,j),refinefactor)*to);
00087 if (!wto.empty()) {
00088 const int& idx = base::U().parentidx(TimeF,Level+1,c,j);
00089 if (_RescaleNonEq)
00090 lbmfixupops::scale_indices(base::U()(TimeC,Level,idx),uh,wto,indices[d],Nindices[d],fac);
00091 else
00092 lbmfixupops::copy_indices(base::U()(TimeC,Level,idx),uh,wto,indices[d],Nindices[d]);
00093 }
00094 }
00095 }
00096 }
00097 end_forall
00098
00099 END_WATCH_FIXUP_WHOLE
00100 START_WATCH
00101 AdaptiveBoundaryUpdate(base::U(),TimeC,Level);
00102 END_WATCH(BOUNDARIES_INTERPOLATION)
00103 START_WATCH
00104 Sync(base::U(),TimeC,Level);
00105 END_WATCH_FIXUP_SYNC
00106 START_WATCH
00107 ExternalBoundaryUpdate(base::U(),TimeC,Level);
00108 END_WATCH(BOUNDARIES_EXTERNAL)
00109 START_WATCH
00110
00111 forall (base::U(),TimeF,Level+1,c)
00112 for (int d=0; d<2*dim; d++) {
00113 if (base::ValidSide(TimeF,Level+1,c,d)) {
00114 BBox corbb(CorrectionBBox(coarsen(base::U().interiorbbox(TimeF,Level+1,c),refinefactor),d));
00115 for (register int j=0; j<base::U().parents(TimeF,Level+1,c); j++) {
00116 const int& idx = base::U().parentidx(TimeF,Level+1,c,j);
00117 BBox to = base::U().interiorbbox(Time,Level,idx)*corbb;
00118 if (!to.empty())
00119 LBM().LocalStep(base::U()(Time,Level,idx),base::U()(TimeC,Level,idx),to,dt);
00120 }
00121 }
00122 }
00123 end_forall
00124 }
00125
00126 BBox CorrectionBBox(const BBox& interior, const int s) {
00127 int d = s/2;
00128 #ifdef DEBUG_PRINT
00129 assert (d>=0 && d<dim);
00130 #endif
00131 BBox bb(interior);
00132 bb.grow(1);
00133 if (s%2 == 0)
00134 bb.upper(d) = bb.lower(d);
00135 else
00136 bb.lower(d) = bb.upper(d);
00137 return bb;
00138 }
00139
00140 inline LBMType& LBM() { return _lbm; }
00141 inline const LBMType& LBM() const { return _lbm; }
00142
00143 protected:
00144 LBMType &_lbm;
00145 int _RescaleNonEq;
00146 int Nindices[2*dim];
00147 int *indices[2*dim];
00148 };
00149
00150 #endif