amroc/amr/ExactSolution.h

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 
00003 // Copyright (C) 2002 Ralf Deiterding
00004 // Brandenburgische Universitaet Cottbus
00005 //
00006 // Copyright (C) 2003-2007 California Institute of Technology
00007 // Ralf Deiterding, ralf@cacr.caltech.edu
00008 
00009 #ifndef AMROC_EXACT_SOLUTION_H
00010 #define AMROC_EXACT_SOLUTION_H
00011 
00025 template <class VectorType, int dim>
00026 class ExactSolution : public AMRBase<VectorType,dim> {
00027   typedef AMRBase<VectorType,dim> base;
00028   typedef typename VectorType::InternalDataType DataType;
00029 
00030 public:
00031   typedef typename base::vec_grid_fct_type vec_grid_fct_type;  
00032   typedef typename base::vec_grid_data_type vec_grid_data_type;
00033   typedef GridFunction<DataType,dim> grid_fct_type;
00034   typedef GridData<DataType,dim> grid_data_type;
00035 
00036   ExactSolution() : base(), ENorm(0), ENormCutOut(1), ENormAll(1), ENormComp(1) {
00037     FileName = "-";
00038   } 
00039 
00040   virtual ~ExactSolution() {}
00041 
00042   virtual void SetGrid(vec_grid_data_type& gd, grid_data_type& gdw, 
00043                        const int& level, const double& t) {}
00044 
00045   virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00046     LocCtrl = Ctrl.getSubDevice(prefix+"Error");    
00047     RegisterAt(LocCtrl,"Norm",ENorm);    
00048     RegisterAt(LocCtrl,"NoRestrict",ENormCutOut);
00049     RegisterAt(LocCtrl,"CombineLevels",ENormAll);
00050     RegisterAt(LocCtrl,"ComponentWise",ENormComp);
00051     RegisterAt(LocCtrl,"FileName",FileName);
00052   }
00053   virtual void register_at(ControlDevice& Ctrl) { register_at(Ctrl, ""); }
00054   virtual void init() {}
00055   virtual void update() {}
00056   virtual void finish() {}
00057 
00058   virtual bool Set(vec_grid_fct_type& u, grid_fct_type& work, const double& t) { 
00059     for (int lev=0; lev<=FineLevel(base::GH()); lev++) {
00060       int Time = CurrentTime(base::GH(),lev); 
00061       int TStep = TimeStep(u,lev);
00062       forall (u,Time+TStep,lev,c) 
00063         u(Time+TStep,lev,c).copy(u(Time,lev,c));
00064         SetGrid(u(Time+TStep,lev,c),work(Time+TStep,lev,c),lev,t);
00065       end_forall      
00066     }
00067     return true;
00068   }    
00069 
00070   //---------------------------------------------------------------------- 
00071   // Calculate Error Norm for conservative variables 
00072   //---------------------------------------------------------------------- 
00073   virtual void ErrorNorm(vec_grid_fct_type& u, grid_fct_type& work, 
00074                          const double& t) {
00075     if (ENorm<1 || ENorm>3) 
00076       return;
00077     if (!Set(u,work,t))
00078       return;
00079 
00080     int me = MY_PROC;
00081     
00082     std::ofstream outfile;
00083     std::ostream* out = 0;
00084     if (me==VizServer && FileName.c_str()[0] != '-') {
00085       outfile.open(FileName.c_str(), std::ios::out | std::ios::app);
00086       out = new std::ostream(outfile.rdbuf());
00087     }
00088     if (me == VizServer) {
00089       std::cout << "                         Error:";
00090       if (out) *out << t;
00091     }
00092     for (int lev=0; lev<=FineLevel(base::GH()); lev++) {
00093       int Time = CurrentTime(base::GH(),lev); 
00094       int TStep = TimeStep(u,lev);
00095       forall (u,Time+TStep,lev,c) 
00096         u(Time+TStep,lev,c).minus(u(Time,lev,c));
00097       end_forall      
00098     }
00099 
00100     DataType *Error = (DataType *) 0;
00101     if (ENormComp) {
00102       for (int i=0; i<base::NEquations(); i++) { 
00103         for (int lev=0; lev<=FineLevel(base::GH()); lev++) {
00104           int Time = CurrentTime(base::GH(),lev); 
00105           int TStep = TimeStep(u,lev);
00106           forall (u,Time+TStep,lev,c) 
00107             equals_from(work(Time+TStep,lev,c), u(Time+TStep,lev,c), i);
00108           end_forall 
00109         }       
00110         CalculateNorm(work, Error);
00111         if (ENormAll) {
00112           std::cout << "   " << Error[0];
00113           if (out) *out << " " << Error[0];
00114         }
00115         else for (int lev=0; lev<=FineLevel(base::GH()); lev++) {
00116           std::cout << "   L(" << i << "," << lev << ")=" << Error[lev];
00117           if (out) *out << " " << Error[lev];     
00118         }
00119       }
00120     }
00121     else {
00122       CalculateNorm(u, Error);
00123       if (ENormAll) {
00124         std::cout << "   " << Error[0];
00125         if (out) *out << " " << Error[0];
00126       }
00127       else for (int lev=0; lev<=FineLevel(base::GH()); lev++) {
00128         std::cout << "   L(" << lev << ")=" << Error[lev];
00129         if (out) *out << " " << Error[lev];       
00130       }
00131     }
00132     if (me == VizServer) {
00133       std::cout << std::endl;
00134       if (out) *out << std::endl;         
00135       outfile.close();
00136       delete out;
00137     }
00138     if (Error) 
00139       delete [] Error;
00140   }  
00141 
00142   virtual void CalculateNorm(vec_grid_fct_type& u, DataType*& Error) {
00143     if (ENormAll)
00144       Error = new DataType[1];
00145     else
00146       Error = new DataType[FineLevel(base::GH())];
00147 
00148     Error[0] = 0.0;
00149     for (int lev=0; lev<=FineLevel(base::GH()); lev++) {
00150       int Time = CurrentTime(base::GH(),lev); 
00151       int TStep = TimeStep(u,lev), FTStep=0;
00152       if (lev<FineLevel(base::GH()))
00153         FTStep = TimeStep(u,lev+1);
00154       
00155       forall (u,Time+TStep,lev,c) 
00156         if (lev<FineLevel(base::GH()) && ENormCutOut) {
00157           forall (u,Time+FTStep,lev+1,c2) 
00158             BBox bb(coarsen(u.interiorbbox(Time+FTStep,lev+1,c2),base::GH().refinefactor(lev)));
00159             u(Time+TStep,lev,c).equals(VectorType(0.0),bb);
00160           end_forall      
00161         }
00162       end_forall      
00163 
00164       DataType ErrorLev = Norm(u,Time+TStep,lev,ENorm);
00165       if (ENormAll) {
00166         switch (ENorm) {
00167         case DAGHNormL1:
00168         case DAGHNormL2:
00169           Error[0] += ErrorLev;
00170           break;
00171         case DAGHNormMax: 
00172           Error[0] = Max(Error[0], ErrorLev);
00173           break;
00174         default:
00175           break;
00176         }
00177       }
00178       else 
00179         Error[lev] = ErrorLev;
00180     }
00181   }  
00182 
00183   virtual void CalculateNorm(grid_fct_type& u, DataType*& Error) {
00184     if (ENormAll)
00185       Error = new DataType[1];
00186     else
00187       Error = new DataType[FineLevel(base::GH())];
00188 
00189     Error[0] = 0.0;
00190     for (int lev=0; lev<=FineLevel(base::GH()); lev++) {
00191       int Time = CurrentTime(base::GH(),lev); 
00192       int TStep = TimeStep(u,lev), FTStep=0;
00193       if (lev<FineLevel(base::GH()))
00194         FTStep = TimeStep(u,lev+1);
00195       
00196       forall (u,Time+TStep,lev,c) 
00197         if (lev<FineLevel(base::GH()) && ENormCutOut) {
00198           forall (u,Time+FTStep,lev+1,c2) 
00199             BBox bb(coarsen(u.interiorbbox(Time+FTStep,lev+1,c2),base::GH().refinefactor(lev)));
00200             u(Time+TStep,lev,c).equals(DataType(0.0),bb);
00201           end_forall      
00202         }
00203       end_forall      
00204 
00205       DataType ErrorLev = Norm(u,Time+TStep,lev,ENorm);
00206       if (ENormAll) {
00207         switch (ENorm) {
00208         case DAGHNormL1:
00209         case DAGHNormL2:
00210           Error[0] += ErrorLev;
00211           break;
00212         case DAGHNormMax: 
00213           Error[0] = Max(Error[0], ErrorLev);
00214           break;
00215         default:
00216           break;
00217         }
00218       }
00219       else 
00220         Error[lev] = ErrorLev;
00221     }
00222   }  
00223 
00224  protected:
00225   int ENorm, ENormCutOut, ENormAll, ENormComp;
00226   std::string FileName;
00227   ControlDevice LocCtrl;  
00228 };
00229 
00230 
00231 #endif
Generated on Thu Jun 30 02:19:17 2016 for AMROC Fluid-solver Framework by 1.6.3