amroc/amr/Criteria/LimiterCriterion.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_LIMITERCRITERION_H
00010 #define AMROC_LIMITERCRITERION_H
00011 
00018 #include "StdCriterion.h"
00019 
00020 #define UROUND 1.0e-15
00021 
00022 template <class VectorType, class FlagType, int dim> class LimiterCriterion; 
00023 
00030 template <class VectorType, class FlagType>
00031 class LimiterCriterion<VectorType,FlagType,1> : 
00032   public StdCriterion<VectorType,FlagType,1> {
00033   typedef typename VectorType::InternalDataType DataType;
00034   typedef StdCriterion<VectorType,FlagType,1> base;
00035 public:
00036   typedef typename base::grid_fct_type grid_fct_type;
00037   typedef typename base::flag_fct_type flag_fct_type;
00038   typedef typename base::grid_data_type grid_data_type;
00039 
00040   LimiterCriterion() : base() {}
00041   
00042   virtual ~LimiterCriterion() {}
00043   
00044   //******************************************************************************
00045   // Abstract class interface
00046   //******************************************************************************
00047   virtual DataType phi(const DataType theta) = 0;
00048   //******************************************************************************
00049 
00050   virtual void Limiter(grid_fct_type& work, const int& Time, const int& Level, 
00051                        int* Offset, const DataType& Sc) {
00052 
00053     int TStep = TimeStep(work,Level);
00054     forall(work, Time, Level, c)
00055       BBox OpBox = work.interiorbbox(Time,Level,c);
00056       Coords& OpBox_stepsize = OpBox.stepsize();
00057       BeginFastIndex1(work, work(Time,Level,c).bbox(), 
00058                       work(Time,Level,c).data(), DataType);
00059       BeginFastIndex1(workout, work(Time+TStep,Level,c).bbox(), 
00060                       work(Time+TStep,Level,c).data(), DataType);
00061 
00062       Coords os = Coords(1,Offset)*OpBox_stepsize;
00063       for_1 (n, OpBox, OpBox_stepsize)
00064         DataType ui  = FastIndex1(work,n);
00065         DataType uim = FastIndex1(work,n-os(0));
00066         DataType uip = FastIndex1(work,n+os(0));
00067         DataType uc = static_cast<DataType>(1.0);
00068         if (std::fabs(ui-uip)>Sc || std::fabs(ui-uim)>Sc) {
00069           if (std::fabs(ui-uip)>UROUND)
00070             uc = phi(std::fabs((ui-uim)/(ui-uip)) );
00071           else uc = static_cast<DataType>(0.0);
00072         }
00073         if (uc < FastIndex1(workout,n))
00074           FastIndex1(workout,n) = uc;   
00075       end_for
00076 
00077       EndFastIndex1(work);
00078       EndFastIndex1(workout);
00079     end_forall
00080   }
00081 };
00082 
00083 
00090 template <class VectorType, class FlagType>
00091 class LimiterCriterion<VectorType,FlagType,2> : 
00092   public StdCriterion<VectorType,FlagType,2> {
00093   typedef typename VectorType::InternalDataType DataType;
00094   typedef StdCriterion<VectorType,FlagType,2> base;
00095 public:
00096   typedef typename base::grid_fct_type grid_fct_type;
00097   typedef typename base::flag_fct_type flag_fct_type;
00098   typedef typename base::grid_data_type grid_data_type;
00099 
00100   LimiterCriterion() : base() {}
00101 
00102   virtual ~LimiterCriterion() {}
00103   
00104   //******************************************************************************
00105   // Abstract class interface
00106   //******************************************************************************
00107   virtual DataType phi(const DataType theta) = 0;
00108   //******************************************************************************
00109 
00110   virtual void Limiter(grid_fct_type& work, const int& Time, const int& Level, 
00111                        int* Offset, const DataType& Sc) {
00112 
00113     int TStep = TimeStep(work,Level);
00114     forall(work, Time, Level, c)
00115       BBox OpBox = work.interiorbbox(Time,Level,c);
00116       Coords& OpBox_stepsize = OpBox.stepsize();
00117       BeginFastIndex2(work, work(Time,Level,c).bbox(), 
00118                       work(Time,Level,c).data(), DataType);
00119       BeginFastIndex2(workout, work(Time+TStep,Level,c).bbox(), 
00120                       work(Time+TStep,Level,c).data(), DataType);
00121 
00122       Coords os = Coords(2,Offset)*OpBox_stepsize;
00123       for_2 (n, m, OpBox, OpBox_stepsize)
00124         DataType ui  = FastIndex2(work,n,m);
00125         DataType uim = FastIndex2(work,n-os(0),m-os(1));
00126         DataType uip = FastIndex2(work,n+os(0),m+os(1));
00127         DataType uc = static_cast<DataType>(1.0);
00128         if (std::fabs(ui-uip)>Sc || std::fabs(ui-uim)>Sc) {
00129           if (std::fabs(ui-uip)>UROUND)
00130             uc = phi(std::fabs((ui-uim)/(ui-uip)) );
00131           else uc = static_cast<DataType>(0.0);
00132         }
00133         if (uc < FastIndex2(workout,n,m))
00134           FastIndex2(workout,n,m) = uc; 
00135       end_for
00136 
00137       EndFastIndex2(work);
00138       EndFastIndex2(workout);
00139     end_forall
00140   }
00141 };
00142 
00143 
00150 template <class VectorType, class FlagType>
00151 class LimiterCriterion<VectorType,FlagType,3> : 
00152   public StdCriterion<VectorType,FlagType,3> {
00153   typedef typename VectorType::InternalDataType DataType;
00154   typedef StdCriterion<VectorType,FlagType,3> base;
00155 public:
00156   typedef typename base::grid_fct_type grid_fct_type;
00157   typedef typename base::flag_fct_type flag_fct_type;
00158   typedef typename base::grid_data_type grid_data_type;
00159 
00160   LimiterCriterion() : base() {}
00161 
00162   virtual ~LimiterCriterion() {}
00163   
00164   //******************************************************************************
00165   // Abstract class interface
00166   //******************************************************************************
00167   virtual DataType phi(const DataType theta) = 0;
00168   //******************************************************************************
00169 
00170   virtual void Limiter(grid_fct_type& work, const int& Time, const int& Level, 
00171                        int* Offset, const DataType& Sc) {
00172 
00173     int TStep = TimeStep(work,Level);
00174     forall(work, Time, Level, c)
00175       BBox OpBox = work.interiorbbox(Time,Level,c);
00176       Coords& OpBox_stepsize = OpBox.stepsize();
00177       BeginFastIndex3(work, work(Time,Level,c).bbox(), 
00178                       work(Time,Level,c).data(), DataType);
00179       BeginFastIndex3(workout, work(Time+TStep,Level,c).bbox(), 
00180                       work(Time+TStep,Level,c).data(), DataType);
00181 
00182       Coords os = Coords(3,Offset)*OpBox_stepsize;
00183       for_3 (n, m, l, OpBox, OpBox_stepsize)
00184         DataType ui  = FastIndex3(work,n,m,l);
00185         DataType uim = FastIndex3(work,n-os(0),m-os(1),l-os(2));
00186         DataType uip = FastIndex3(work,n+os(0),m+os(1),l+os(2));
00187         DataType uc = static_cast<DataType>(1.0);
00188         if (std::fabs(ui-uip)>Sc || std::fabs(ui-uim)>Sc) {
00189           if (std::fabs(ui-uip)>UROUND)
00190             uc = phi(std::fabs((ui-uim)/(ui-uip)) );
00191           else uc = static_cast<DataType>(0.0);
00192         }
00193         if (uc < FastIndex3(workout,n,m,l))
00194           FastIndex3(workout,n,m,l) = uc;       
00195       end_for
00196 
00197       EndFastIndex3(work);
00198       EndFastIndex3(workout);
00199     end_forall
00200   }
00201 };
00202 
00203 
00204 #endif
Generated on Thu Jun 30 02:19:14 2016 for AMROC Fluid-solver Framework by 1.6.3