amroc/clawpack/ClpFixup.h

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 
00003 // Copyright (C) 2002 Ralf Deiterding
00004 // Brandenburgische Universitaet Cottbus
00005 
00006 #ifndef AMROC_CLP_FIXUP_H
00007 #define AMROC_CLP_FIXUP_H
00008 
00016 #include "AMRFixup.h"
00017 
00018 // -*- C++ -*-
00019 
00020 #include <iostream>
00021 #include <cstdlib>
00022 
00023 template <class VectorType, class AuxVectorType, int dim> class ClpIntegrator;
00024 
00031 template <class VectorType, int dim>
00032 class ClpFixupVecOps : private FixupOps<VectorType,VectorType,dim> {
00033 public:
00034   ClpFixupVecOps() {}
00035 protected:
00036   inline void copyv_to(GridData<VectorType,minus_1<dim>::dim> &target, 
00037                        const GridData<VectorType,dim> &source, 
00038                        const BBox &where, const int s) 
00039     { FixupOps<VectorType,VectorType,dim>::copy_to(target, source, where, s); }
00040 };
00041 
00042 #define f_rcflx FORTRAN_NAME(rcflx, RCFLX)
00043 extern "C" {
00044   void f_rcflx(const INTEGER& mflx);
00045 }
00046 
00053 template <class VectorType, class FixupType, class AuxVectorType, int dim>
00054 class ClpFixup : 
00055   public AMRFixup<VectorType,FixupType,dim>,
00056   public ClpFixupVecOps<VectorType,dim> {
00057 
00058   typedef typename VectorType::InternalDataType DataType;
00059   typedef AMRFixup<VectorType,FixupType,dim> base;
00060   typedef ClpFixupVecOps<VectorType,dim> ops;
00061 
00062 public:
00063   typedef ClpIntegrator<VectorType,AuxVectorType,dim> integrator_type; 
00064   typedef typename base::vec_grid_data_type vec_grid_data_type;
00065   typedef typename base::ld_vec_grid_data_type ld_vec_grid_data_type;
00066 
00067   ClpFixup(integrator_type& integ) : base(), _Integrator(integ) {}
00068 
00069 protected:
00070   void CorrectFirstOrder(const int Time, const int Level, const int c, 
00071                          const double tc, const double tf,
00072                          const double dtc, const double dtf, 
00073                          const int& mdim) {
00074     
00075     int TimeC = CurrentTime(base::GH(),Level-1);
00076     DCoords dx = base::GH().worldStep(base::U().GF_StepSize(Level));
00077     const int refinefactor = RefineFactor(base::GH(),Level-1);
00078 #ifdef DEBUG_PRINT_WHOLE_FIXUP
00079     ( comm_service::log() << "\n*** Adding First Order Fluxes - "
00080       << base::U().interiorbbox(Time,Level,c) << " ***\n").flush();
00081 #endif
00082       
00083     double Fac[2*dim]; 
00084     int d;
00085     for (d=0; d<dim; d++) {
00086       Fac[2*d] = dtf; Fac[2*d+1] = -dtf;
00087       for (int dd=0; dd<dim; dd++)
00088         if (d != dd) {
00089           Fac[2*d] *= dx(dd); Fac[2*d+1] *= dx(dd);
00090         }
00091     }
00092 
00093     for (d=0; d<2*dim; d++) 
00094       if ((mdim==0 || d/2==mdim-1) && base::ValidSide(Time,Level,c,d)) {
00095 #ifdef DEBUG_PRINT_WHOLE_FIXUP
00096         ( comm_service::log() << "Side:  " << d << "\n").flush();
00097 #endif    
00098         BBox bb(base::StateVecBBox(base::U().interiorbbox(Time,Level,c), d));
00099         BBox bbC(base::StateVecBBox(coarsen(base::U().interiorbbox(Time,Level,c),
00100                                       refinefactor),d));      
00101         BBox workbb(bb); base::AlignBBox(workbb, d);
00102         
00103         ld_vec_grid_data_type LeftState(workbb); 
00104         ld_vec_grid_data_type RightState(workbb); 
00105         ld_vec_grid_data_type NewState(workbb); 
00106         BBox bbl(bb), bbr(bb);
00107         double dtl = dtf, dtr = dtf, 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(base::U()(Time,Level,c), Level, bb, tf,
00115                "ClpFixup2::CorrectFirstOrder::LeftState,fine");
00116 #endif
00117           ops::copyv_to(LeftState, base::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<base::U().parents(Time,Level,c) ; j++) {
00122             BBox wbbC(coarsen(base::U().parentbox(Time,Level,c,j),refinefactor)*bbC);
00123             if (!wbbC.empty()) {
00124 #ifdef DEBUG_AMR
00125               if (Integrator_().Abort()) 
00126                 Integrator_().CheckGrid(base::U()(TimeC,Level-1,base::U().parentidx(Time,Level,c,j)), Level-1, 
00127                                         wbbC, tc, "ClpFixup2::CorrectFirstOrder::RightState,coarse");
00128 #endif
00129               ops::copyv_to(RightState, base::U()(TimeC,Level-1,base::U().parentidx(Time,Level,c,j)), wbbC, d);
00130             }
00131           }
00132           tr = tc;
00133           dtr = dtc; 
00134           bbr.coarsen(d/2, refinefactor);
00135         }
00136         else {
00137           /* Coarse grid data comes from all parent grids */
00138           for (register int j=0; j<base::U().parents(Time,Level,c) ; j++) {
00139             BBox wbbC(coarsen(base::U().parentbox(Time,Level,c,j),refinefactor)*bbC);
00140             if (!wbbC.empty()) {
00141 #ifdef DEBUG_AMR
00142               if (Integrator_().Abort()) 
00143                 Integrator_().CheckGrid(base::U()(TimeC,Level-1,base::U().parentidx(Time,Level,c,j)), Level-1, 
00144                                         wbbC, tc, "ClpFixup2::CorrectFirstOrder::LeftState,coarse");
00145 #endif
00146               ops::copyv_to(LeftState, base::U()(TimeC,Level-1,base::U().parentidx(Time,Level,c,j)), wbbC, d);
00147             }
00148           }
00149           tl = tc;
00150           dtl = dtc; 
00151           bbl.coarsen(d/2, refinefactor);
00152         
00153           /* Fine grid data comes from the current grid */
00154 #ifdef DEBUG_AMR
00155           if (Integrator_().Abort()) 
00156             Integrator_().CheckGrid(base::U()(Time,Level,c), Level, bb, tf,
00157                "ClpFixup2::CorrectFirstOrder::RightState,fine");
00158 #endif
00159           ops::copyv_to(RightState, base::U()(Time,Level,c), bb, d);
00160           tr = tf;
00161         }
00162       
00163         ((integrator_type&)Integrator_()).CalculateRP(LeftState, bbl, tl, dtl,  
00164            RightState, bbr, tr, dtr, NewState, d/2+1);      
00165       
00166 #ifdef DEBUG_PRINT_WHOLE_FIXUP
00167         base::debug_print_ld(base::F(d)(Time,Level,c));
00168         ( comm_service::log() << "\n").flush();
00169         ( comm_service::log() << "LeftState: \n").flush();
00170         base::debug_print_ldv(LeftState); ( comm_service::log() << "\n").flush();
00171         ( comm_service::log() << "RightState: \n").flush();
00172         base::debug_print_ldv(RightState); ( comm_service::log() << "\n").flush();
00173         ( comm_service::log() << "NewState: \n").flush();
00174         base::debug_print_ldv(NewState); ( comm_service::log() << "\n").flush();
00175 #endif
00176         NewState.multiply(Fac[d],workbb);
00177         base::addf_to(base::F(d)(Time,Level,c),NewState,workbb);
00178 #ifdef DEBUG_PRINT_WHOLE_FIXUP
00179         base::debug_print_ld(base::F(d)(Time,Level,c));
00180         ( comm_service::log() << "\n").flush();
00181 #endif  
00182       }
00183   } 
00184   
00185   virtual void AddFluxes(const int Time, const int Level, const int c,
00186                          vec_grid_data_type* flux[], 
00187                          const double tc, const double tf, 
00188                          const double dtc, const double dtf, 
00189                          const int& mdim) {
00190     int mflx;
00191     f_rcflx(mflx);
00192     if (mflx==0) {
00193       if (dim>1 && Integrator_().DimSplit()) {
00194         std::cerr << "\n\nThe fixup is incorrect, if dimensional splitting is used with a\n";
00195         std::cerr << "scheme returning flux differences. Set Method(3)>=0 instead.\n";
00196         std::cerr << "Aborting programm...\n" << std::flush;
00197         std::exit(-1);
00198       }
00199       CorrectFirstOrder(Time, Level, c, tc, tf, dtc, dtf, mdim);
00200     }
00201     
00202     base::AddIntercellFluxes(Time, Level, c, flux, dtf, mdim);
00203   }
00204 
00205   inline integrator_type& Integrator_() { return _Integrator; }
00206   inline const integrator_type& Integrator_() const { return _Integrator; }
00207   inline const int& NGhosts() const { return _Integrator.NGhosts(); }
00208 
00209 private:
00210   integrator_type& _Integrator;
00211 }; 
00212 
00213 #endif