AMROC - Source code documentation

00001 #ifndef AMROC_AMRFIXUP_H
00002 #define AMROC_AMRFIXUP_H
00003 
00011 #if DIM == 1
00012 # include "Fixup1.h"
00013 #elif DIM == 2
00014 # include "Fixup2.h"
00015 #elif DIM == 3  
00016 # include "Fixup3.h"
00017 #endif
00018 
00019 #include "AMRFixupInterface.h"
00020 
00021 
00022 #ifndef AMRFixupName
00023 #define AMRFixup(dim)      name2(AMRFixup,dim)
00024 #define AMRFixupName
00025 #endif
00026 
00047 template <class VectorType, class FixupType>
00048 class AMRFixup(DIM) : public AMRFixupInterface(DIM)<VectorType,FixupType>,
00049     public Fixup(DIM)<VectorType,FixupType>,
00050     private FixupOps(DIM)<VectorType,FixupType> {
00051 
00052   typedef Integrator(DIM)<VectorType> integrator_type;    
00053   typedef typename VectorType::InternalDataType DataType;
00054   typedef GridData(DIM)<VectorType>  vec_grid_data_type;
00055 
00056   typedef GridData(DIM_1)<VectorType> ld_vec_grid_data_type;
00057   typedef GridData(DIM_1)<FixupType>  ld_fixup_grid_data_type;
00058 
00059   typedef GridFunction(DIM)<VectorType> vec_grid_fct_type;  
00060   typedef GridFunction(DIM_1)<FixupType> ld_fixup_grid_fct_type;  
00061 
00062 public:
00063   AMRFixup(DIM)(integrator_type& integ) : 
00064     Fixup(DIM)<VectorType,FixupType>(integ) {}
00065 
00066   virtual void register_at(ControlDevice& Ctrl, const string& prefix) {
00067     LocCtrl = Ctrl.getSubDevice(prefix+"AMRFixup");
00068     RegisterAt(LocCtrl,"Equations",_FixupEquations);
00069   } 
00070   virtual void register_at(ControlDevice& Ctrl) { register_at(Ctrl, ""); }
00071 
00072   virtual void SaveFluxes(const int Time, const int Level, const int c,
00073                           vec_grid_data_type* flux[], 
00074                           const double dt, DCoords& dx, 
00075                           const int& mdim) {
00076 
00077     if (!U().has_children(Time,Level,c)) return;
00078     const int refinefactor = RefineFactor(GH(),Level);
00079     DataType Fac = -dt;
00080     int d;
00081     for (d=0; d<DIM; d++) 
00082       Fac *= dx(d);
00083     
00084 #ifdef DEBUG_PRINT_WHOLE_FIXUP
00085     ( comm_service::log() << "\n*** Saving Coarse Fluxes - " 
00086       << U().interiorbbox(Time,Level,c) << " ***\n").flush();
00087 #endif
00088     for (register int j=0; j<U().children(Time,Level,c); j++) {
00089       BBox work(U().childbox(Time,Level,c,j)*U().interiorbbox(Time,Level,c));
00090       if (!work.empty()) {
00091         const int& idx = U().childidx(Time,Level,c,j);
00092         for (d=0; d<2*DIM; d++) 
00093           if ((mdim==0 || d/2==mdim-1) && ValidSide(Time,Level+1,idx,d)) {
00094             BBox bb(FixupBBox(coarsen(U().interiorbbox(Time,Level+1,idx),refinefactor),d)*work);
00095             if (!bb.empty()) {
00096               bb = FluxBBox(bb,d);
00097 #ifdef DEBUG_PRINT_WHOLE_FIXUP
00098               ( comm_service::log() << "To " << U().interiorbbox(Time,Level+1,idx) 
00099                 << " Side: " << d << " using " << bb << " into " 
00100                 << F(d)(Time,Level+1,idx).bbox() << "\n").flush();
00101 #endif
00102               (*flux[d]).multiply(Fac/dx(d/2), bb);
00103               copyf_to(F(d)(Time,Level+1,idx),(*flux[d]),bb,d);
00104               (*flux[d]).divide(Fac/dx(d/2), bb);
00105 #ifdef DEBUG_PRINT_WHOLE_FIXUP
00106               AlignBBox(bb,d);
00107               debug_print_ld(F(d)(Time,Level+1,idx),bb); 
00108               ( comm_service::log() << "\n" ).flush();
00109 #endif
00110             }
00111           }
00112       }
00113     }
00114   }
00115 
00116   virtual void AddFluxes(const int Time, const int Level, const int c,
00117                          vec_grid_data_type* flux[], 
00118                          const double tc, const double tf, 
00119                          const double dt, DCoords& dx, 
00120                          const int& mdim) {
00121     int mflx;
00122     f_rcflx(mflx);
00123     if (mflx==0) {
00124       cerr << "AMRCLAWFixup must be used instead of AMRFixup!\nAborting programm...\n" << flush;
00125 #ifdef DEBUG_PRINT
00126       ( comm_service::log() << "AMRCLAWFixup must be used instead of AMRFixup!"
00127         << "\nAborting programm...\n\n" ).flush();
00128 #endif
00129       assert(0);
00130     }
00131     AddIntercellFluxes(Time, Level, c, flux, dt, dx, mdim);
00132   }
00133   
00134   virtual void Correction(const int Time, const int WTime, const int Level, DCoords& dx) {
00135     int t_sten = 0;   
00136     int s_sten = 1;       
00137     char GFName[DAGHBktGFNameWidth];
00138     sprintf(GFName,"u-fixup");
00139     vec_grid_fct_type uh(GFName, t_sten, s_sten, CurrentTime(GH(),Level), Level, 
00140                          GH(), DAGHCellCentered, 0, 1, DAGHCommSimple, 
00141                          DAGHNoBoundary, DAGHNoAdaptBoundary, DAGHNoExternalGhost);
00142 
00143 #ifdef DEBUG_PRINT_WHOLE_FIXUP
00144     ( comm_service::log() << "\n*********** Conservative Fixup ***********\n").flush();
00145 #endif
00146 
00147     for (int d=0; d<2*DIM; d++) 
00148       ParallelFixup(Time, Level, dx, d, uh);
00149   }  
00150 
00151 protected:  
00152   int ValidSide(const int Time, const int Level, const int c, const int s) const {
00153     return (U().has_adaptiveboundary(Time,Level,c,s));
00154   }
00155     
00156   void AddIntercellFluxes(const int Time, const int Level, const int c,
00157                           vec_grid_data_type* flux[], const double dt,
00158                           DCoords& dx, const int& mdim) {
00159     
00160     DataType Fac = dt;
00161     int d;
00162     for (d=0; d<DIM; d++) 
00163       Fac *= dx(d);
00164 
00165 #ifdef DEBUG_PRINT_WHOLE_FIXUP
00166     ( comm_service::log() << "\n*** Adding Higher Order Fluxes - " 
00167       << U().interiorbbox(Time,Level,c) << " ***\n").flush();
00168 #endif
00169     for (d=0; d<2*DIM; d++) 
00170       if ((mdim==0 || d/2==mdim-1) && ValidSide(Time,Level,c,d)) {
00171         BBox bb(FluxBBox(U().interiorbbox(Time,Level,c), d));
00172         (*flux[d]).multiply(Fac/dx(d/2), bb);
00173         addf_to(F(d)(Time,Level,c),(*flux[d]),bb,d);
00174         (*flux[d]).divide(Fac/dx(d/2), bb);
00175 #ifdef DEBUG_PRINT_WHOLE_FIXUP
00176         debug_print_ld(F(d)(Time,Level,c));
00177         ( comm_service::log() << "\n").flush();
00178 #endif
00179       }
00180   }
00181   
00182   void ParallelFixup(const int Time, const int Level, 
00183                      DCoords& dx, const int d, vec_grid_fct_type& uh) {
00184 
00185     int TimeF = CurrentTime(GH(),Level+1);
00186     int TimeC = CurrentTime(GH(),Level);
00187     const int refinefactor = RefineFactor(GH(),Level);
00188     DataType Fac = 1.0;
00189     for (int dd=0; dd<DIM; dd++) 
00190       Fac *= dx(dd);
00191 
00192     forall (uh,TimeC,Level,c)  
00193       uh(TimeC,Level,c) = 0.0;
00194     end_forall
00195 
00196     forall (U(),TimeF,Level+1,c) 
00197       if (ValidSide(TimeF,Level+1,c,d)) {
00198         BBox to(FixupBBox(coarsen(U().interiorbbox(TimeF,Level+1,c),refinefactor),d));
00199         for (register int j=0; j<U().parents(TimeF,Level+1,c) ; j++) {
00200           BBox wto(coarsen(U().parentbox(TimeF,Level+1,c,j),refinefactor)*to);
00201           const int& idx = U().parentidx(TimeF,Level+1,c,j);
00202           if (!wto.empty()) {
00203             copyf_from(uh(TimeC,Level,idx),wto,F(d)(TimeF,Level+1,c),d);
00204             uh(TimeC,Level,idx).multiply((d%2==0 ? -1.0 : 1.0) / Fac,wto);
00205 #ifdef DEBUG_PRINT_WHOLE_FIXUP
00206             ( comm_service::log() << "Writing at Side: " << d << " of " 
00207               << U().interiorbbox(TimeF,Level+1,c) << " using " << wto << "\n").flush();
00208             debug_print(uh(TimeC,Level,idx), wto);
00209             ( comm_service::log() << "\n" ).flush(); 
00210 #endif  
00211           }
00212         }
00213       }
00214     end_forall 
00215 
00216     END_WATCH_FIXUP_WHOLE
00217     START_WATCH
00218       Sync(uh,TimeC,Level);
00219     END_WATCH_FIXUP_SYNC
00220     START_WATCH
00221 
00222     forall (U(),Time,Level,c)  
00223       BBox to = U().interiorbbox(Time,Level,c);
00224       BBox from(to);
00225       from.shift(d/2, d%2 ? -1 : 1);
00226       U()(Time,Level,c).plus(uh(TimeC,Level,c), to, from);
00227     end_forall 
00228   }
00229 
00230   BBox StateVecBBox(const BBox& interior, const int s) {    
00231     int d = s/2;
00232 #ifdef DEBUG_PRINT
00233     assert (d>=0 && d<DIM);
00234 #endif
00235     BBox bb(interior); 
00236     if (s%2 == 0) {
00237       bb.growlower(d,-1);
00238       bb.upper(d) = bb.lower(d);
00239     }
00240     else {
00241       bb.growupper(d,1);
00242       bb.lower(d) = bb.upper(d);
00243     }
00244     return bb;
00245   }
00246 
00247   BBox FixupBBox(const BBox& interior, const int s) {    
00248     int d = s/2;
00249 #ifdef DEBUG_PRINT
00250     assert (d>=0 && d<DIM);
00251 #endif
00252     BBox bb(interior); 
00253     if (s%2 == 0) 
00254       bb.upper(d) = bb.lower(d);
00255     else 
00256       bb.lower(d) = bb.upper(d);
00257     return bb;
00258   }
00259 
00260   virtual BBox FluxBBox(const BBox& interior, const int s) {    
00261     int d = s/2;
00262 #ifdef DEBUG_PRINT
00263     assert (d>=0 && d<DIM);
00264 #endif
00265     BBox bb(interior); 
00266     if (s%2 == 0)
00267       bb.upper(d) = bb.lower(d);
00268     else {
00269       bb.growupper(d,1);
00270       bb.lower(d) = bb.upper(d);
00271     }
00272     return bb;
00273   }
00274 
00275 protected:
00276   inline void copyf_to(ld_fixup_grid_data_type &target, const BBox &to, 
00277                        const vec_grid_data_type &source, 
00278                        const BBox &from, const int s) 
00279     { copy_to(target, to, source, from, s); }
00280   inline void copyf_to(ld_fixup_grid_data_type &target, 
00281                        const vec_grid_data_type &source, 
00282                        const BBox &where, const int s) 
00283     { copy_to(target, source, where, s); }
00284 
00285   inline void copyf_from(vec_grid_data_type &target, const BBox &to, 
00286                          const ld_fixup_grid_data_type &source, 
00287                          const BBox &from, const int s) 
00288     { copy_from(target, to, source, from, s); }
00289   inline void copyf_from(vec_grid_data_type &target, const BBox &where, 
00290                          const ld_fixup_grid_data_type &source, const int s) 
00291     { copy_from(target, where, source, s); }
00292 
00293   inline void addf_to(ld_fixup_grid_data_type &target, const BBox &to, 
00294                       const vec_grid_data_type &source, 
00295                       const BBox &from, const int s) 
00296     { add_to(target, to, source, from, s); }
00297   inline void addf_to(ld_fixup_grid_data_type &target, 
00298                       const vec_grid_data_type &source, 
00299                       const BBox &where, const int s) 
00300     { add_to(target, source, where, s); }
00301 
00302   inline void addf_from(vec_grid_data_type &target, const BBox &to, 
00303                          const ld_fixup_grid_data_type &source, 
00304                          const BBox &from, const int s) 
00305     { add_from(target, to, source, from, s); }
00306   inline void addf_from(vec_grid_data_type &target, const BBox &where, 
00307                          const ld_fixup_grid_data_type &source, const int s) 
00308     { add_from(target, where, source, s); }
00309 
00310   inline void addf_to(ld_fixup_grid_data_type &target, const BBox &to, 
00311                        const ld_vec_grid_data_type &source, 
00312                        const BBox &from) 
00313     { add_to(target, to, source, from); }  
00314   inline void addf_to(ld_fixup_grid_data_type &target, 
00315                        const ld_vec_grid_data_type &source, 
00316                        const BBox &where) 
00317     { add_to(target, source, where); }  
00318 
00319 };
00320 
00321 #endif
AMRFixup.h

last update: 06/01/04
