AMROC - Source code documentation

00001 #ifndef AMROC_CLP_FIXUP_H
00002 #define AMROC_CLP_FIXUP_H
00003 
00011 #if DIM == 1
00012 # include "ClpIntegrator1.h"
00013 #elif DIM == 2
00014 # include "ClpIntegrator2.h"
00015 #elif DIM == 3  
00016 # include "ClpIntegrator3.h"
00017 #endif
00018 #include "AMRFixup.h"
00019 
00020 #ifndef ClpFixupVecOpsName
00021 #define ClpFixupVecOps(dim)      name2(ClpFixupVecOps,dim)
00022 #define ClpFixupVecOpsName
00023 #endif
00024 
00031 template <class VectorType>
00032 class ClpFixupVecOps(DIM) : private FixupOps(DIM)<VectorType,VectorType> {
00033 public:
00034   ClpFixupVecOps(DIM)() {}
00035 protected:
00036   inline void copyv_to(GridData(DIM_1)<VectorType> &target, 
00037                        const GridData(DIM)<VectorType> &source, 
00038                        const BBox &where, const int s) 
00039     { copy_to(target, source, where, s); }
00040 };
00041 
00042 #ifndef ClpFixupName
00043 #define ClpFixup(dim)      name2(ClpFixup,dim)
00044 #define ClpFixupName
00045 #endif
00046 
00053 template <class VectorType, class FixupType, class AuxVectorType>
00054 class ClpFixup(DIM) : 
00055   public AMRFixup(DIM)<VectorType,FixupType>,
00056   public ClpFixupVecOps(DIM)<VectorType> {
00057 
00058   typedef ClpIntegrator(DIM)<VectorType,AuxVectorType> integrator_type; 
00059   typedef AMRFixup(DIM)<VectorType,FixupType> fixup_base_type;
00060 
00061   typedef typename VectorType::InternalDataType DataType;
00062   typedef GridData(DIM)<VectorType>  vec_grid_data_type;
00063   typedef GridData(DIM_1)<VectorType> ld_vec_grid_data_type;
00064 
00065 public:
00066   ClpFixup(DIM)(integrator_type& integ) : fixup_base_type(integ) {}
00067 
00068 protected:
00069   void CorrectFirstOrder(const int Time, const int Level, const int c, 
00070                          const double tc, const double tf,
00071                          const double dt, DCoords& dx, 
00072                          const int& mdim) {
00073     
00074     int TimeC = CurrentTime(GH(),Level-1);
00075     const int refinefactor = RefineFactor(GH(),Level-1);
00076 #ifdef DEBUG_PRINT_WHOLE_FIXUP
00077     ( comm_service::log() << "\n*** Adding First Order Fluxes - "
00078       << U().interiorbbox(Time,Level,c) << " ***\n").flush();
00079 #endif
00080       
00081     double Fac[2*DIM]; 
00082     int d;
00083     for (d=0; d<DIM; d++) {
00084       Fac[2*d] = dt; Fac[2*d+1] = -dt;
00085       for (int dd=0; dd<DIM; dd++)
00086         if (d != dd) {
00087           Fac[2*d] *= dx(dd); Fac[2*d+1] *= dx(dd);
00088         }
00089     }
00090 
00091     for (d=0; d<2*DIM; d++) 
00092       if ((mdim==0 || d/2==mdim-1) && ValidSide(Time,Level,c,d)) {
00093 #ifdef DEBUG_PRINT_WHOLE_FIXUP
00094         ( comm_service::log() << "Side:  " << d << "\n").flush();
00095 #endif    
00096         BBox bb(StateVecBBox(U().interiorbbox(Time,Level,c), d));
00097         BBox bbC(StateVecBBox(coarsen(U().interiorbbox(Time,Level,c),
00098                                       refinefactor),d));      
00099         BBox workbb(bb); AlignBBox(workbb, d);
00100         
00101         ld_vec_grid_data_type LeftState(workbb); 
00102         ld_vec_grid_data_type RightState(workbb); 
00103         ld_vec_grid_data_type NewState(workbb); 
00104         DCoords dxl(dx); 
00105         DCoords dxr(dx); 
00106         BBox bbl(bb), bbr(bb);
00107         double dtl = dt, dtr = dt, tl, tr;
00108       
00109         /* Collect data for Riemann problems at coarse-fine boundaries */
00110         if (d % 2 == 0) {
00111           /* Fine grid data comes from the current grid */
00112 #ifdef DEBUG_AMR
00113           if (Integrator_().Abort()) 
00114             Integrator_().CheckGrid(U()(Time,Level,c), bb, 
00115                "ClpFixup2::CorrectFirstOrder::LeftState,fine");
00116 #endif
00117           copyv_to(LeftState, U()(Time,Level,c), bb, d);
00118           tl = tf;
00119 
00120           /* Coarse grid data comes from all parent grids */
00121           for (register int j=0; j<U().parents(Time,Level,c) ; j++) {
00122             BBox wbbC(coarsen(U().parentbox(Time,Level,c,j),refinefactor)*bbC);
00123             if (!wbbC.empty()) {
00124 #ifdef DEBUG_AMR
00125               if (Integrator_().Abort()) 
00126                 Integrator_().CheckGrid(U()(TimeC,Level-1,U().parentidx(Time,Level,c,j)), wbbC, 
00127                    "ClpFixup2::CorrectFirstOrder::RightState,coarse");
00128 #endif
00129               copyv_to(RightState, U()(TimeC,Level-1,U().parentidx(Time,Level,c,j)), wbbC, d);
00130             }
00131           }
00132           tr = tc;
00133           dtr      *= refinefactor; 
00134           dxr(d/2) *= refinefactor;
00135           bbr.coarsen(d/2, refinefactor);
00136         }
00137         else {
00138           /* Coarse grid data comes from all parent grids */
00139           for (register int j=0; j<U().parents(Time,Level,c) ; j++) {
00140             BBox wbbC(coarsen(U().parentbox(Time,Level,c,j),refinefactor)*bbC);
00141             if (!wbbC.empty()) {
00142 #ifdef DEBUG_AMR
00143               if (Integrator_().Abort()) 
00144                 Integrator_().CheckGrid(U()(TimeC,Level-1,U().parentidx(Time,Level,c,j)), wbbC, 
00145                    "ClpFixup2::CorrectFirstOrder::LeftState,coarse");
00146 #endif
00147               copyv_to(LeftState, U()(TimeC,Level-1,U().parentidx(Time,Level,c,j)), wbbC, d);
00148             }
00149           }
00150           tl = tc;
00151           dtl      *= refinefactor; 
00152           dxl(d/2) *= refinefactor;         
00153           bbl.coarsen(d/2, refinefactor);
00154         
00155           /* Fine grid data comes from the current grid */
00156 #ifdef DEBUG_AMR
00157           if (Integrator_().Abort()) 
00158             Integrator_().CheckGrid(U()(Time,Level,c), bb, 
00159                "ClpFixup2::CorrectFirstOrder::RightState,fine");
00160 #endif
00161           copyv_to(RightState, U()(Time,Level,c), bb, d);
00162           tr = tf;
00163         }
00164       
00165         ((integrator_type&)Integrator_()).CalculateRP(LeftState, bbl, tl, dtl, dxl, 
00166            RightState, bbr, tr, dtr, dxr, NewState, d/2+1);         
00167       
00168 #ifdef DEBUG_PRINT_WHOLE_FIXUP
00169         debug_print_ld(F(d)(Time,Level,c));
00170         ( comm_service::log() << "\n").flush();
00171         ( comm_service::log() << "LeftState: \n").flush();
00172         debug_print_ldv(LeftState); ( comm_service::log() << "\n").flush();
00173         ( comm_service::log() << "RightState: \n").flush();
00174         debug_print_ldv(RightState); ( comm_service::log() << "\n").flush();
00175         ( comm_service::log() << "NewState: \n").flush();
00176         debug_print_ldv(NewState); ( comm_service::log() << "\n").flush();
00177 #endif
00178         NewState.multiply(Fac[d],workbb);
00179         addf_to(F(d)(Time,Level,c),NewState,workbb);
00180 #ifdef DEBUG_PRINT_WHOLE_FIXUP
00181         debug_print_ld(F(d)(Time,Level,c));
00182         ( comm_service::log() << "\n").flush();
00183 #endif  
00184       }
00185   } 
00186   
00187   virtual void AddFluxes(const int Time, const int Level, const int c,
00188                          vec_grid_data_type* flux[], 
00189                          const double tc, const double tf, 
00190                          const double dt, DCoords& dx,
00191                          const int& mdim) {
00192     int mflx;
00193     f_rcflx(mflx);
00194     if (mflx==0)
00195       CorrectFirstOrder(Time, Level, c, tc, tf, dt, dx, mdim);
00196     
00197     AddIntercellFluxes(Time, Level, c, flux, dt, dx, mdim);
00198   }
00199 }; 
00200 
00201 #endif
ClpFixup.h

last update: 06/01/04
