File List
amroc/converters/Evaluators/SCompEulerEvaluator.h

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_SCOMP_EULER_EVALUATOR_H
00010 #define AMROC_SCOMP_EULER_EVALUATOR_H
00011 
00019 #include <iostream>
00020 #include <string>
00021 #include <cstdio>
00022 #include <cmath>
00023 #include <cfloat>
00024 
00025 #include "Evaluator.h"
00026 
00033 template <class VectorType>
00034 class SCompEulerEvaluator : public Evaluator<VectorType> {
00035   typedef Evaluator<VectorType> base;
00036   typedef typename VectorType::InternalDataType DataType;
00037 
00038 public:
00039   typedef typename base::vector_block_type vector_block_type;
00040   typedef typename base::data_block_type data_block_type;
00041 
00042   SCompEulerEvaluator() : base(), RU(1.), PA(1.), Gamma(1.), W(1.) {
00043     std::sprintf(base::title,"3D-Euler equations - Single ideal gas");
00044   }
00045 
00046   virtual void update() {   
00047     ControlDevice ChemCtrl(GetFileControlDevice("chem.dat",""));
00048     RegisterAt(ChemCtrl,"RU",RU); 
00049     RegisterAt(ChemCtrl,"PA",PA); 
00050     RegisterAt(ChemCtrl,"Gamma",Gamma); 
00051     RegisterAt(ChemCtrl,"Sp",SpName); 
00052     RegisterAt(ChemCtrl,"W",W); 
00053     ChemCtrl.update();  
00054   }
00055 
00056   virtual void SetUp(VisualGridBase *gr) {
00057     base::SetUp(gr);
00058     const int bk=base::NKeys();
00059     std::sprintf(&(base::tkeys[bk    *LEN_TKEYS]), "Temperature             ");
00060     std::sprintf(&(base::tkeys[(bk+1)*LEN_TKEYS]), "Pressure                ");
00061     std::sprintf(&(base::tkeys[(bk+2)*LEN_TKEYS]), "Entropy                 ");
00062     std::sprintf(&(base::tkeys[(bk+3)*LEN_TKEYS]), "Speed of Sound          ");
00063     std::sprintf(&(base::tkeys[(bk+4)*LEN_TKEYS]), "Machnumber              ");
00064     std::sprintf(&(base::tkeys[(bk+5)*LEN_TKEYS]), "Schl.-Plot Temperature  ");
00065     std::sprintf(&(base::tkeys[(bk+6)*LEN_TKEYS]), "Schl.-Plot Pressure     ");
00066     std::sprintf(&(base::tkeys[(bk+7)*LEN_TKEYS]), "Schl.-Plot Entropy      ");
00067     std::sprintf(&(base::tkeys[(bk+8)*LEN_TKEYS]), "Schl.-Plot Sound-Speed  ");
00068     std::sprintf(&(base::tkeys[(bk+9)*LEN_TKEYS]), "Schl.-Plot Machnumber   ");
00069     base::fkeys[bk  ] = base::Grid->FKeys();
00070     base::fkeys[bk+1] = base::Grid->FKeys();
00071     base::fkeys[bk+2] = base::Grid->FKeys();
00072     base::fkeys[bk+3] = base::Grid->FKeys();
00073     base::fkeys[bk+4] = base::Grid->FKeys();
00074     base::fkeys[bk+5] = base::Grid->FKeys();
00075     base::fkeys[bk+6] = base::Grid->FKeys();
00076     base::fkeys[bk+7] = base::Grid->FKeys();
00077     base::fkeys[bk+8] = base::Grid->FKeys();
00078     base::fkeys[bk+9] = base::Grid->FKeys();
00079     base::ikeys[bk  ] = 116; base::ikeys[bk+1] = 112; base::ikeys[bk+2] = 110; 
00080     base::ikeys[bk+3] = 99;  base::ikeys[bk+4] = 109; 
00081     base::ikeys[bk+5] = 84;  base::ikeys[bk+6] = 80;  base::ikeys[bk+7] = 78;  
00082     base::ikeys[bk+8] = 67;  base::ikeys[bk+9] = 77;  
00083   }
00084 
00085   virtual int NKeys() const { return base::NKeys()+10; }
00086 
00087   virtual void SetScalCells(int key,float v[],int& offset, 
00088                             vector_block_type& DB, float* dx, bool* CompRead) {
00089     int idx=0;
00090     int keym1 = key-1;
00091     const int bk=base::NKeys();
00092     BeginFastIndex3(dat, DB.bbox(), DB.data(), VectorType);
00093     if (key>=bk+1 && key<=bk+10) {        
00094       if (!CompRead[4] || !CompRead[0])
00095         return;
00096 
00097       float rho, u1, u2, u3, E, p;
00098       for_3 (i,j,k,DB.bbox(),DB.bbox().stepsize())
00099         rho = FastIndex3(dat,i,j,k)(0);               
00100         u1  = FastIndex3(dat,i,j,k)(1);
00101         u2  = FastIndex3(dat,i,j,k)(2);
00102         u3  = FastIndex3(dat,i,j,k)(3);
00103         E   = FastIndex3(dat,i,j,k)(4);  
00104         p   = (Gamma-1.0)*(E-0.5*rho*(u1*u1+u2*u2+u3*u3));
00105         if      (key==bk+1 || key==bk+6)
00106           v[idx+offset] = (p*W)/(rho*RU);
00107         else if (key==bk+2 || key==bk+7)
00108           v[idx+offset] = p;
00109         else if (key==bk+3 || key==bk+8) 
00110           v[idx+offset] = RU/(W*(Gamma-1))*std::log(p/std::pow(rho,Gamma));
00111         else if (key==bk+4 || key==bk+9)
00112           v[idx+offset] = std::sqrt(Gamma*p/rho);
00113         else if (key==bk+5 || key==bk+10) 
00114           v[idx+offset] = std::sqrt(u1*u1+u2*u2+u3*u3)/std::sqrt(Gamma*p/rho);
00115 
00116         base::mvals[2*keym1]   = Min(base::mvals[2*keym1],   v[idx+offset]);
00117         base::mvals[2*keym1+1] = Max(base::mvals[2*keym1+1], v[idx+offset]);
00118         idx++;        
00119       end_for
00120     }
00121     else if (key<=bk)
00122        base::SetScalCells(key,v,offset,DB,dx,CompRead);
00123     
00124     EndFastIndex3(dat); 
00125     offset+=idx;
00126 
00127     if (key>=bk+6 && key<=bk+10)
00128       std::cerr << &(base::tkeys[keym1*LEN_TKEYS])
00129                 << " is NOT supported for Type=6. Use Type=1 instead!"
00130                 << std::endl; 
00131   }
00132     
00133   virtual void SetScalNodes(int key,float v[],int& offset,
00134                             vector_block_type& DB, const BBox& bbox, 
00135                             float* dx, bool* CompRead) {
00136     int idx=0;
00137     int keym1 = key-1;
00138     const int bk=base::NKeys();
00139     BeginFastIndex3(dat, DB.bbox(), DB.data(), VectorType);
00140     const Coords& step = bbox.stepsize();
00141     if (key>=bk+1 && key<=bk+5) {         
00142       if (!CompRead[4] || !CompRead[0])
00143         return;
00144 
00145       float rho, u1, u2, u3, E, p, c;
00146       for_3 (i, j, k, bbox, step)
00147         BBox AddOver(3,i-step(0),j-step(1),k-step(2),
00148                      i,j,k,step(0),step(1),step(2));
00149         AddOver *= DB.bbox();
00150           
00151         c=0;
00152         v[idx+offset] = 0.0;
00153         for_3 (l, m, n, AddOver, step)
00154           rho = FastIndex3(dat,l,m,n)(0);                     
00155           u1  = FastIndex3(dat,l,m,n)(1);
00156           u2  = FastIndex3(dat,l,m,n)(2);
00157           u3  = FastIndex3(dat,l,m,n)(3);
00158           E   = FastIndex3(dat,l,m,n)(4);
00159           if (rho<FLT_MAX && u1<FLT_MAX && u2<FLT_MAX && 
00160               u3<FLT_MAX && E<FLT_MAX) { 
00161             p = (Gamma-1.0)*(E-0.5*rho*(u1*u1+u2*u2+u3*u3));
00162             if      (key==bk+1)
00163               v[idx+offset] += (p*W)/(rho*RU);
00164             else if (key==bk+2)
00165               v[idx+offset] += p;
00166             else if (key==bk+3)
00167               v[idx+offset] += RU/(W*(Gamma-1))*
00168                                std::log(p/std::pow(rho,Gamma));
00169             else if (key==bk+4)
00170               v[idx+offset] += std::sqrt(Gamma*p/rho);
00171             else if (key==bk+5)
00172               v[idx+offset] += std::sqrt(u1*u1+u2*u2+u3*u3)/
00173                                std::sqrt(Gamma*p/rho);
00174             c += 1.0;
00175           }
00176         end_for
00177         if (c>0) v[idx+offset] /= c;
00178         else v[idx+offset] = FLT_MAX;
00179 
00180         base::mvals[2*keym1]   = Min(base::mvals[2*keym1],   v[idx+offset]);
00181         base::mvals[2*keym1+1] = Max(base::mvals[2*keym1+1], v[idx+offset]);
00182         idx++;        
00183       end_for 
00184     }
00185 
00186     if (key>=bk+6 && key<=bk+10) {        
00187       if (!CompRead[4] || !CompRead[0])
00188         return;
00189 
00190       float rho, u1, u2, u3, E, p;
00191       for_3 (i, j, k, bbox, step)
00192         BBox AddOver(3,i-step(0),j-step(1),k-step(2),
00193                      i,j,k,step(0),step(1),step(2));
00194         AddOver *= DB.bbox();
00195           
00196         data_block_type DBHelp(AddOver);
00197         DBHelp = FLT_MAX;
00198         BeginFastIndex3(help, DBHelp.bbox(), DBHelp.data(), DataType);
00199         for_3 (l, m, n, AddOver, step)
00200           rho = FastIndex3(dat,l,m,n)(0);                     
00201           u1  = FastIndex3(dat,l,m,n)(1);
00202           u2  = FastIndex3(dat,l,m,n)(2);
00203           u3  = FastIndex3(dat,l,m,n)(3);
00204           E   = FastIndex3(dat,l,m,n)(4);
00205 
00206           if (rho<FLT_MAX && u1<FLT_MAX && u2<FLT_MAX && 
00207               u3<FLT_MAX && E<FLT_MAX) {
00208             p = (Gamma-1.0)*(E-0.5*rho*(u1*u1+u2*u2+u3*u3));
00209             if      (key==bk+6)
00210               FastIndex3(help,l,m,n) = (p*W)/(rho*RU);
00211             else if (key==bk+7)
00212               FastIndex3(help,l,m,n) = p;
00213             else if (key==bk+8)
00214               FastIndex3(help,l,m,n) = RU/(W*(Gamma-1))*
00215                                        std::log(p/std::pow(rho,Gamma));
00216             else if (key==bk+9)
00217               FastIndex3(help,l,m,n) = std::sqrt(Gamma*p/rho);
00218             else if (key==bk+10)
00219               FastIndex3(help,l,m,n) = std::sqrt(u1*u1+u2*u2+u3*u3)/
00220                                        std::sqrt(Gamma*p/rho);
00221           }
00222         end_for
00223         EndFastIndex3(help); 
00224         v[idx+offset] = base::SetScalGradNodes(DBHelp,dx);
00225 
00226         base::mvals[2*keym1]   = Min(base::mvals[2*keym1],   v[idx+offset]);
00227         base::mvals[2*keym1+1] = Max(base::mvals[2*keym1+1], v[idx+offset]);
00228         idx++;        
00229       end_for 
00230     }
00231 
00232     if (key<=bk) 
00233       base::SetScalNodes(key,v,offset,DB,bbox,dx,CompRead);
00234 
00235     EndFastIndex3(dat); 
00236     offset+=idx;
00237   }
00238 
00239 protected:
00240   DataType RU, PA, Gamma, W; 
00241   std::string SpName;
00242 };
00243 
00244 #endif
Generated on Thu Jun 30 02:19:02 2016 for AMROC's HDF Tools by 1.6.3