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