vtf-logo

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@amroc.net
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         if (uc < FastIndex1(workout,n))
00073           FastIndex1(workout,n) = uc;   
00074       end_for
00075 
00076       EndFastIndex1(work);
00077       EndFastIndex1(workout);
00078     end_forall
00079   }
00080 };
00081 
00082 
00089 template <class VectorType, class FlagType>
00090 class LimiterCriterion<VectorType,FlagType,2> : 
00091   public StdCriterion<VectorType,FlagType,2> {
00092   typedef typename VectorType::InternalDataType DataType;
00093   typedef StdCriterion<VectorType,FlagType,2> base;
00094 public:
00095   typedef typename base::grid_fct_type grid_fct_type;
00096   typedef typename base::flag_fct_type flag_fct_type;
00097   typedef typename base::grid_data_type grid_data_type;
00098 
00099   LimiterCriterion() : base() {}
00100 
00101   virtual ~LimiterCriterion() {}
00102   
00103   //******************************************************************************
00104   // Abstract class interface
00105   //******************************************************************************
00106   virtual DataType phi(const DataType theta) = 0;
00107   //******************************************************************************
00108 
00109   virtual void Limiter(grid_fct_type& work, const int& Time, const int& Level, 
00110                        int* Offset, const DataType& Sc) {
00111 
00112     int TStep = TimeStep(work,Level);
00113     forall(work, Time, Level, c)
00114       BBox OpBox = work.interiorbbox(Time,Level,c);
00115       Coords& OpBox_stepsize = OpBox.stepsize();
00116       BeginFastIndex2(work, work(Time,Level,c).bbox(), 
00117                       work(Time,Level,c).data(), DataType);
00118       BeginFastIndex2(workout, work(Time+TStep,Level,c).bbox(), 
00119                       work(Time+TStep,Level,c).data(), DataType);
00120 
00121       Coords os = Coords(2,Offset)*OpBox_stepsize;
00122       for_2 (n, m, OpBox, OpBox_stepsize)
00123         DataType ui  = FastIndex2(work,n,m);
00124         DataType uim = FastIndex2(work,n-os(0),m-os(1));
00125         DataType uip = FastIndex2(work,n+os(0),m+os(1));
00126         DataType uc = static_cast<DataType>(1.0);
00127         if (std::fabs(ui-uip)>Sc || std::fabs(ui-uim)>Sc)
00128           if (std::fabs(ui-uip)>UROUND)
00129             uc = phi(std::fabs((ui-uim)/(ui-uip)) );
00130           else uc = static_cast<DataType>(0.0);
00131         if (uc < FastIndex2(workout,n,m))
00132           FastIndex2(workout,n,m) = uc; 
00133       end_for
00134 
00135       EndFastIndex2(work);
00136       EndFastIndex2(workout);
00137     end_forall
00138   }
00139 };
00140 
00141 
00148 template <class VectorType, class FlagType>
00149 class LimiterCriterion<VectorType,FlagType,3> : 
00150   public StdCriterion<VectorType,FlagType,3> {
00151   typedef typename VectorType::InternalDataType DataType;
00152   typedef StdCriterion<VectorType,FlagType,3> base;
00153 public:
00154   typedef typename base::grid_fct_type grid_fct_type;
00155   typedef typename base::flag_fct_type flag_fct_type;
00156   typedef typename base::grid_data_type grid_data_type;
00157 
00158   LimiterCriterion() : base() {}
00159 
00160   virtual ~LimiterCriterion() {}
00161   
00162   //******************************************************************************
00163   // Abstract class interface
00164   //******************************************************************************
00165   virtual DataType phi(const DataType theta) = 0;
00166   //******************************************************************************
00167 
00168   virtual void Limiter(grid_fct_type& work, const int& Time, const int& Level, 
00169                        int* Offset, const DataType& Sc) {
00170 
00171     int TStep = TimeStep(work,Level);
00172     forall(work, Time, Level, c)
00173       BBox OpBox = work.interiorbbox(Time,Level,c);
00174       Coords& OpBox_stepsize = OpBox.stepsize();
00175       BeginFastIndex3(work, work(Time,Level,c).bbox(), 
00176                       work(Time,Level,c).data(), DataType);
00177       BeginFastIndex3(workout, work(Time+TStep,Level,c).bbox(), 
00178                       work(Time+TStep,Level,c).data(), DataType);
00179 
00180       Coords os = Coords(3,Offset)*OpBox_stepsize;
00181       for_3 (n, m, l, OpBox, OpBox_stepsize)
00182         DataType ui  = FastIndex3(work,n,m,l);
00183         DataType uim = FastIndex3(work,n-os(0),m-os(1),l-os(2));
00184         DataType uip = FastIndex3(work,n+os(0),m+os(1),l+os(2));
00185         DataType uc = static_cast<DataType>(1.0);
00186         if (std::fabs(ui-uip)>Sc || std::fabs(ui-uim)>Sc)
00187           if (std::fabs(ui-uip)>UROUND)
00188             uc = phi(std::fabs((ui-uim)/(ui-uip)) );
00189           else uc = static_cast<DataType>(0.0);
00190         if (uc < FastIndex3(workout,n,m,l))
00191           FastIndex3(workout,n,m,l) = uc;       
00192       end_for
00193 
00194       EndFastIndex3(work);
00195       EndFastIndex3(workout);
00196     end_forall
00197   }
00198 };
00199 
00200 
00201 #endif

Generated on Fri Aug 24 13:00:52 2007 for AMROC Fluid-solver Framework - by  doxygen 1.4.7