00001
00002
00003
00004
00005
00006
00007
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
00179 base::mvals[2*keym1] = Min(base::mvals[2*keym1], v[idx+offset]);
00180 base::mvals[2*keym1+1] = Max(base::mvals[2*keym1+1], v[idx+offset]);
00181 idx++;
00182 end_for
00183 }
00184
00185 if (key>=bk+6 && key<=bk+10) {
00186 if (!CompRead[4] || !CompRead[0])
00187 return;
00188
00189 float rho, u1, u2, u3, E, p;
00190 for_3 (i, j, k, bbox, step)
00191 BBox AddOver(3,i-step(0),j-step(1),k-step(2),
00192 i,j,k,step(0),step(1),step(2));
00193 AddOver *= DB.bbox();
00194
00195 data_block_type DBHelp(AddOver);
00196 DBHelp = FLT_MAX;
00197 BeginFastIndex3(help, DBHelp.bbox(), DBHelp.data(), DataType);
00198 for_3 (l, m, n, AddOver, step)
00199 rho = FastIndex3(dat,l,m,n)(0);
00200 u1 = FastIndex3(dat,l,m,n)(1);
00201 u2 = FastIndex3(dat,l,m,n)(2);
00202 u3 = FastIndex3(dat,l,m,n)(3);
00203 E = FastIndex3(dat,l,m,n)(4);
00204
00205 if (rho<FLT_MAX && u1<FLT_MAX && u2<FLT_MAX &&
00206 u3<FLT_MAX && E<FLT_MAX) {
00207 p = (Gamma-1.0)*(E-0.5*rho*(u1*u1+u2*u2+u3*u3));
00208 if (key==bk+6)
00209 FastIndex3(help,l,m,n) = (p*W)/(rho*RU);
00210 else if (key==bk+7)
00211 FastIndex3(help,l,m,n) = p;
00212 else if (key==bk+8)
00213 FastIndex3(help,l,m,n) = RU/(W*(Gamma-1))*
00214 std::log(p/std::pow(rho,Gamma));
00215 else if (key==bk+9)
00216 FastIndex3(help,l,m,n) = std::sqrt(Gamma*p/rho);
00217 else if (key==bk+10)
00218 FastIndex3(help,l,m,n) = std::sqrt(u1*u1+u2*u2+u3*u3)/
00219 std::sqrt(Gamma*p/rho);
00220 }
00221 end_for
00222 EndFastIndex3(help);
00223 v[idx+offset] = base::SetScalGradNodes(DBHelp,dx);
00224
00225 base::mvals[2*keym1] = Min(base::mvals[2*keym1], v[idx+offset]);
00226 base::mvals[2*keym1+1] = Max(base::mvals[2*keym1+1], v[idx+offset]);
00227 idx++;
00228 end_for
00229 }
00230
00231 if (key<=bk)
00232 base::SetScalNodes(key,v,offset,DB,bbox,dx,CompRead);
00233
00234 EndFastIndex3(dat);
00235 offset+=idx;
00236 }
00237
00238 protected:
00239 DataType RU, PA, Gamma, W;
00240 std::string SpName;
00241 };
00242
00243 #endif
