00001
00002
00003
00004
00005
00006 #ifndef AMROC_PROBLEM_H
00007 #define AMROC_PROBLEM_H
00008
00009 #include "eulerrhok3.h"
00010
00011 #define NEQUATIONS 15
00012 #define NEQUSED 14
00013 #define NFIXUP 13
00014 #define NAUX 6
00015
00016 #include "ClpProblem.h"
00017
00018 #define OWN_FLAGGING
00019 #define OWN_AMRSOLVER
00020 #include "ClpStdProblem.h"
00021 #include "F77Interfaces/F77UpdateLevelTransfer.h"
00022 #include "SpecialOutput/AMRGridAccumulation.h"
00023 #include "AMRInterpolation.h"
00024
00025 #define f_rcveloc FORTRAN_NAME(rcveloc, RCVELOC)
00026 #define f_track FORTRAN_NAME(track, TRACK)
00027 extern "C" {
00028 void f_rcveloc(const DOUBLE v[]);
00029 void f_track(const INTEGER &, const DOUBLE *, const DOUBLE *, DOUBLE *, INTEGER&,
00030 const DOUBLE&, const DOUBLE&, const DOUBLE&);
00031 }
00032
00033 class FlaggingSpecific :
00034 public AMRFlagging<VectorType,FixupType,FlagType,DIM> {
00035 typedef AMRFlagging<VectorType,FixupType,FlagType,DIM> base;
00036 public:
00037 FlaggingSpecific(base::solver_type& solver) : base(solver) {
00038 base::AddCriterion(new ScaledGradient<VectorType,FlagType,DIM>());
00039 base::AddCriterion(new LimiterType<VectorType,FlagType,DIM>());
00040 base::AddCriterion(new AbsoluteError<VectorType,FixupType,FlagType,DIM>(solver));
00041 base::AddCriterion(new RelativeError<VectorType,FixupType,FlagType,DIM>(solver));
00042 base::AddCriterion(new F77ScaledGradient<VectorType,FlagType,DIM>(f_flg));
00043 base::AddCriterion(new F77LimiterType<VectorType,FlagType,DIM>(f_flg));
00044 base::AddCriterion(new F77AbsoluteError<VectorType,FixupType,FlagType,DIM>(solver,f_flg));
00045 base::AddCriterion(new F77RelativeError<VectorType,FixupType,FlagType,DIM>(solver,f_flg));
00046 }
00047 ~FlaggingSpecific() { DeleteAllCriterions(); }
00048
00049 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00050 base::register_at(Ctrl, prefix);
00051 char VariableName[32];
00052 for (int d=0; d<DIM; d++) {
00053 sprintf(VariableName,"DCellsMin(%d)",d+1);
00054 RegisterAt(base::LocCtrl,VariableName,dcmin[d]);
00055 sprintf(VariableName,"DCellsMax(%d)",d+1);
00056 RegisterAt(base::LocCtrl,VariableName,dcmax[d]);
00057 }
00058 RegisterAt(base::LocCtrl,"DCellsVeloc",dcv);
00059 }
00060 virtual void register_at(ControlDevice& Ctrl) {
00061 register_at(Ctrl, "");
00062 }
00063
00064 virtual void SetFlags(const int Time, const int Level, double t, double dt) {
00065 base::SetFlags(Time, Level, t, dt);
00066 double dx = DeltaX(base::GH(),0,0);
00067 int dcmaxuse[DIM];
00068 for (int d=0; d<base::Dim(); d++) dcmaxuse[d] = dcmax[d];
00069 dcmaxuse[0] += int(dcv*t/dx);
00070 Coords lb(base::Dim(), dcmin), ub(base::Dim(), dcmaxuse), s(base::Dim(),1);
00071 BBox bbox(lb*RefinedBy(base::GH(),MaxLevel(base::GH())),
00072 ub*RefinedBy(base::GH(),MaxLevel(base::GH())),
00073 s*StepSize(base::GH(),Level));
00074 forall (base::Flags(),Time,Level,c)
00075 base::Flags()(Time,Level,c).equals(GoodPoint, bbox);
00076 end_forall
00077 }
00078
00079 private:
00080 int dcmin[DIM], dcmax[DIM];
00081 double dcv;
00082 };
00083
00084
00085 class SolverSpecific :
00086 public AMRSolver<VectorType,FixupType,FlagType,DIM> {
00087 typedef VectorType::InternalDataType DataType;
00088 typedef AMRSolver<VectorType,FixupType,FlagType,DIM> base;
00089 typedef AMRGridAccumulation<DataType,DIM> accumulation_type;
00090 typedef F77FileOutput<VectorType,DIM> output_type;
00091 typedef AMRInterpolation<VectorType,DIM> interpolation_type;
00092 typedef interpolation_type::point_type point_type;
00093 public:
00094 SolverSpecific(IntegratorSpecific& integ,
00095 base::initial_condition_type& init,
00096 base::boundary_conditions_type& bc) : base(integ, init, bc) {
00097 SetLevelTransfer(new F77UpdateLevelTransfer<VectorType,DIM>(f_prolong, f_restrict, f_tupdate));
00098 SetFileOutput(new F77FileOutput<VectorType,DIM>(f_out));
00099 SetFixup(new FixupSpecific(integ));
00100 SetFlagging(new FlaggingSpecific(*this));
00101 _Accumulation = new accumulation_type(f_rcveloc);
00102 _Interpolation = new interpolation_type();
00103 MaxType = 1;
00104 TrackFront = 0;
00105 TrackYPos = 0.;
00106 TrackZPos = 0.;
00107 TrackMin = 0.1;
00108 TrackMax = 1.e4;
00109 TrackdfMin = 1000.;
00110 }
00111
00112 ~SolverSpecific() {
00113 delete _LevelTransfer;
00114 delete _Flagging;
00115 delete _Fixup;
00116 delete _FileOutput;
00117 delete _Accumulation;
00118 delete _Interpolation;
00119 }
00120
00121 virtual void register_at(ControlDevice& Ctrl) { base::register_at(Ctrl); }
00122 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00123 base::register_at(Ctrl, prefix);
00124 RegisterAt(base::LocCtrl,"MaxType",MaxType);
00125 RegisterAt(base::LocCtrl,"TrackFront",TrackFront);
00126 RegisterAt(base::LocCtrl,"TrackYPos",TrackYPos);
00127 RegisterAt(base::LocCtrl,"TrackZPos",TrackZPos);
00128 RegisterAt(base::LocCtrl,"TrackMin",TrackMin);
00129 RegisterAt(base::LocCtrl,"TrackMax",TrackMax);
00130 RegisterAt(base::LocCtrl,"TrackdfMin",TrackdfMin);
00131 _Accumulation->register_at(base::LocCtrl, prefix);
00132 }
00133
00134 virtual void SetupData() {
00135 base::SetupData();
00136 _Accumulation->SetupData(base::PGH(), base::NGhosts());
00137 _Interpolation->SetupData(base::PGH(), base::NGhosts());
00138 }
00139
00140 virtual void Initialize_(const double& dt_start) {
00141 base::Initialize_(dt_start);
00142 _Accumulation->Initialize();
00143 }
00144
00145 virtual void Output() {
00146 int Time = CurrentTime(base::GH(),0);
00147 if (Time == base::LastOutputTime()) return;
00148 _Accumulation->GridCombine(VizServer);
00149 int me = MY_PROC;
00150 if (me == VizServer) {
00151 int Level = _Accumulation->AccLevel();
00152 base::FileOutput_().WriteOut(_Accumulation->AccGrid(),_Accumulation->AccName(),
00153 Time,Level,base::t[0]);
00154 }
00155 base::Output();
00156 }
00157
00158 virtual void Checkpointing_(const char* CheckpointFile) {
00159 base::Checkpointing_(CheckpointFile);
00160 _Accumulation->GridCombine(VizServer);
00161 if (MY_PROC == VizServer)
00162 _Accumulation->Checkpointing();
00163 }
00164
00165 virtual bool Restart_(const char* CheckpointFile) {
00166 if (base::Restart_(CheckpointFile)) {
00167 if (MY_PROC == VizServer)
00168 _Accumulation->Restart();
00169 else
00170 _Accumulation->Initialize();
00171 return true;
00172 }
00173 else
00174 return false;
00175 }
00176
00177 virtual void AdvanceLevel(const int Level, int RegridEvery, bool RegridDone,
00178 bool ShadowAllowed, bool DoFixup,
00179 bool RecomposeBaseLev, bool RecomposeHighLev) {
00180
00181 base::AdvanceLevel(Level,RegridEvery,RegridDone,ShadowAllowed,
00182 DoFixup,RecomposeBaseLev,RecomposeHighLev);
00183
00184 int Time = CurrentTime(base::GH(),Level);
00185 if (MaxType) {
00186 AdaptiveBoundaryUpdate(base::U(),Time,Level);
00187 Sync(base::U(),Time,Level);
00188 ExternalBoundaryUpdate(base::U(),Time,Level);
00189
00190 int irho = 1;
00191 int iu = NEQUSED-DIM-2;
00192 int iv = NEQUSED-DIM-1;
00193 int iw = NEQUSED-DIM-0;
00194 int TStep = TimeStep(base::U(),Level);
00195 ((output_type*) base::_FileOutput)->Transform(base::U(), base::Work(), Time+TStep,
00196 Level, irho, base::t[Level]);
00197 forall (base::U(),Time,Level,c)
00198 Coords ss(base::U()(Time,Level,c).bbox().stepsize());
00199 BBox bbi(base::U().interiorbbox(Time,Level,c));
00200 BeginFastIndex3(U, base::U()(Time,Level,c).bbox(),
00201 base::U()(Time,Level,c).data(), VectorType);
00202 BeginFastIndex3(Work, Work()(Time,Level,c).bbox(),
00203 Work()(Time,Level,c).data(), DataType);
00204 BeginFastIndex3(WorkN, Work()(Time+TStep,Level,c).bbox(),
00205 Work()(Time+TStep,Level,c).data(), DataType);
00206 DataType uy, uz, vx, vz, wx, wy, vorx, vory, vorz;
00207 DCoords dx = base::GH().worldStep(ss);
00208 for_3 (i, j, k, bbi, ss)
00209 uy = (FastIndex3(U,i,j+ss(1),k)(iu)/FastIndex3(WorkN,i,j+ss(1),k)-
00210 FastIndex3(U,i,j-ss(1),k)(iu)/FastIndex3(WorkN,i,j-ss(1),k))/(2.0*dx(1));
00211 uz = (FastIndex3(U,i,j,k+ss(2))(iu)/FastIndex3(WorkN,i,j,k+ss(2))-
00212 FastIndex3(U,i,j,k-ss(2))(iu)/FastIndex3(WorkN,i,j,k-ss(2)))/(2.0*dx(2));
00213 vx = (FastIndex3(U,i+ss(0),j,k)(iv)/FastIndex3(WorkN,i+ss(0),j,k)-
00214 FastIndex3(U,i-ss(0),j,k)(iv)/FastIndex3(WorkN,i-ss(0),j,k))/(2.0*dx(0));
00215 vz = (FastIndex3(U,i,j,k+ss(2))(iv)/FastIndex3(WorkN,i,j,k+ss(2))-
00216 FastIndex3(U,i,j,k-ss(2))(iv)/FastIndex3(WorkN,i,j,k-ss(2)))/(2.0*dx(2));
00217 wx = (FastIndex3(U,i+ss(0),j,k)(iw)/FastIndex3(WorkN,i+ss(0),j,k)-
00218 FastIndex3(U,i-ss(0),j,k)(iw)/FastIndex3(WorkN,i-ss(0),j,k))/(2.0*dx(0));
00219 wy = (FastIndex3(U,i,j+ss(1),k)(iw)/FastIndex3(WorkN,i,j+ss(1),k)-
00220 FastIndex3(U,i,j-ss(1),k)(iw)/FastIndex3(WorkN,i,j-ss(1),k))/(2.0*dx(1));
00221 vorx = wy-vz; vory = uz-wx; vorz = vx-uy;
00222 FastIndex3(Work,i,j,k) = std::sqrt(vorx*vorx+vory*vory+vorz*vorz);
00223 end_for
00224 EndFastIndex3(U);
00225 EndFastIndex3(Work);
00226 EndFastIndex3(WorkN);
00227 end_forall
00228 }
00229 else {
00230 int ipress = base::Dim()+4;
00231 ((output_type*) base::_FileOutput)->Transform(base::U(), base::Work(), Time, Level,
00232 ipress, base::t[Level]);
00233 }
00234
00235 _Accumulation->Accumulate(base::Work(), Time, Level, base::t[Level]);
00236 }
00237
00238 virtual double Tick(int VariableTimeStepping, const double dtv[],
00239 const double cflv[], int& Rejections) {
00240
00241 double cfl = base::Tick(VariableTimeStepping,dtv,cflv,Rejections);
00242 if (!TrackFront) return cfl;
00243
00244 int me = MY_PROC;
00245 int Npoints, finelevel = FineLevel(base::GH());
00246 Npoints = base::shape[0]*RefinedBy(base::GH(),finelevel);
00247 double dx = DeltaX(base::GH(),0,finelevel);
00248
00249 DataType* p = new DataType[Npoints];
00250 DataType* r = new DataType[Npoints];
00251 DataType* xr = new DataType[Npoints];
00252 point_type* xc = new point_type[Npoints];
00253 register int n;
00254
00255 for (n=0; n<Npoints; n++) {
00256 r[n] = base::geom[0]+(n+0.5)*dx;
00257 xc[n](0) = r[n]; xc[n](1) = TrackYPos; xc[n](2) = TrackZPos;
00258 }
00259
00260 for (register int l=0; l<=FineLevel(base::GH()); l++) {
00261 int Time = CurrentTime(base::GH(),l);
00262 int press_cnt = base::Dim()+4;
00263 ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time, l,
00264 press_cnt, base::t[l]);
00265 }
00266
00267 _Interpolation->PointsValues(base::Work(),base::t[0],Npoints,xc,p);
00268 _Interpolation->ArrayCombine(VizServer,Npoints,p);
00269
00270 if (me == VizServer) {
00271 int nc;
00272 f_track(Npoints,p,r,xr,nc,TrackMin,TrackMax,TrackdfMin);
00273 std::ofstream outfile;
00274 std::ostream* out;
00275 std::string fname = "xpos.dat";
00276 outfile.open(fname.c_str(), std::ios::out | std::ios::app);
00277 out = new std::ostream(outfile.rdbuf());
00278 *out << base::t[0] << " " << xr[nc-1] << std::endl;
00279 outfile.close();
00280 delete out;
00281 }
00282
00283 delete [] p;
00284 delete [] r;
00285 delete [] xr;
00286 delete [] xc;
00287
00288 return cfl;
00289 }
00290
00291 protected:
00292 int MaxType, TrackFront;
00293 double TrackYPos, TrackZPos, TrackMin, TrackMax, TrackdfMin;
00294 accumulation_type* _Accumulation;
00295 interpolation_type* _Interpolation;
00296 };