00001
00002
00003
00004
00005
00006 #ifndef AMROC_PROBLEM_H
00007 #define AMROC_PROBLEM_H
00008
00009 #include "eulerrhok2.h"
00010
00011 #define NEQUATIONS 14
00012 #define NEQUSED 13
00013 #define NFIXUP 12
00014 #define NAUX 4
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(in2eurhok, IN2EURHOK)
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
00061 class SolverSpecific :
00062 public AMRSolver<VectorType,FixupType,FlagType,DIM> {
00063 typedef VectorType::InternalDataType DataType;
00064 typedef AMRSolver<VectorType,FixupType,FlagType,DIM> base;
00065 typedef AMRGridAccumulation<DataType,DIM> accumulation_type;
00066 typedef F77FileOutput<VectorType,DIM> output_type;
00067 typedef AMRInterpolation<VectorType,DIM> interpolation_type;
00068 typedef interpolation_type::point_type point_type;
00069 public:
00070 SolverSpecific(IntegratorSpecific& integ,
00071 base::initial_condition_type& init,
00072 base::boundary_conditions_type& bc) : base(integ, init, bc) {
00073 SetLevelTransfer(new F77UpdateLevelTransfer<VectorType,DIM>(f_prolong, f_restrict, f_tupdate));
00074 SetFileOutput(new F77FileOutput<VectorType,DIM>(f_out));
00075 SetFixup(new FixupSpecific(integ));
00076 SetFlagging(new FlaggingSpecific(*this));
00077 _Accumulation = new accumulation_type(f_rcveloc);
00078 _Interpolation = new interpolation_type();
00079 MaxType = 1;
00080 TrackFront = 0;
00081 TrackYPos = 0.;
00082 TrackMin = 0.1;
00083 TrackMax = 1.e4;
00084 TrackdfMin = 1000.;
00085 }
00086
00087 ~SolverSpecific() {
00088 delete _LevelTransfer;
00089 delete _Flagging;
00090 delete _Fixup;
00091 delete _FileOutput;
00092 delete _Accumulation;
00093 delete _Interpolation;
00094 }
00095
00096 virtual void register_at(ControlDevice& Ctrl) { base::register_at(Ctrl); }
00097 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00098 base::register_at(Ctrl, prefix);
00099 RegisterAt(base::LocCtrl,"MaxType",MaxType);
00100 RegisterAt(base::LocCtrl,"TrackFront",TrackFront);
00101 RegisterAt(base::LocCtrl,"TrackYPos",TrackYPos);
00102 RegisterAt(base::LocCtrl,"TrackMin",TrackMin);
00103 RegisterAt(base::LocCtrl,"TrackMax",TrackMax);
00104 RegisterAt(base::LocCtrl,"TrackdfMin",TrackdfMin);
00105 _Accumulation->register_at(base::LocCtrl, prefix);
00106 }
00107
00108 virtual void SetupData() {
00109 base::SetupData();
00110 _Accumulation->SetupData(base::PGH(), base::NGhosts());
00111 _Interpolation->SetupData(base::PGH(), base::NGhosts());
00112 }
00113
00114 virtual void Initialize_(const double& dt_start) {
00115 base::Initialize_(dt_start);
00116 _Accumulation->Initialize();
00117 }
00118
00119 virtual void Output() {
00120 int Time = CurrentTime(base::GH(),0);
00121 if (Time == base::LastOutputTime()) return;
00122 _Accumulation->GridCombine(VizServer);
00123 int me = MY_PROC;
00124 if (me == VizServer) {
00125 int Level = _Accumulation->AccLevel();
00126 base::FileOutput_().WriteOut(_Accumulation->AccGrid(),_Accumulation->AccName(),
00127 Time,Level,base::t[0]);
00128 }
00129 base::Output();
00130 }
00131
00132 virtual void Checkpointing_(const char* CheckpointFile) {
00133 base::Checkpointing_(CheckpointFile);
00134 _Accumulation->GridCombine(VizServer);
00135 if (MY_PROC == VizServer)
00136 _Accumulation->Checkpointing();
00137 }
00138
00139 virtual bool Restart_(const char* CheckpointFile) {
00140 if (base::Restart_(CheckpointFile)) {
00141 if (MY_PROC == VizServer)
00142 _Accumulation->Restart();
00143 else
00144 _Accumulation->Initialize();
00145 return true;
00146 }
00147 else
00148 return false;
00149 }
00150
00151 virtual void AdvanceLevel(const int Level, int RegridEvery, bool RegridDone,
00152 bool ShadowAllowed, bool DoFixup,
00153 bool RecomposeBaseLev, bool RecomposeHighLev) {
00154
00155 base::AdvanceLevel(Level,RegridEvery,RegridDone,ShadowAllowed,
00156 DoFixup,RecomposeBaseLev,RecomposeHighLev);
00157
00158 int Time = CurrentTime(base::GH(),Level);
00159 if (MaxType) {
00160 AdaptiveBoundaryUpdate(base::U(),Time,Level);
00161 Sync(base::U(),Time,Level);
00162 ExternalBoundaryUpdate(base::U(),Time,Level);
00163
00164 int irho = 1;
00165 int iu = NEQUSED-DIM-2;
00166 int iv = NEQUSED-DIM-1;
00167 int TStep = TimeStep(base::U(),Level);
00168 ((output_type*) base::_FileOutput)->Transform(base::U(), base::Work(), Time+TStep,
00169 Level, irho, base::t[Level]);
00170 forall (base::U(),Time,Level,c)
00171 Coords ss(base::U()(Time,Level,c).bbox().stepsize());
00172 BBox bbi(base::U().interiorbbox(Time,Level,c));
00173 BeginFastIndex2(U, base::U()(Time,Level,c).bbox(),
00174 base::U()(Time,Level,c).data(), VectorType);
00175 BeginFastIndex2(Work, Work()(Time,Level,c).bbox(),
00176 Work()(Time,Level,c).data(), DataType);
00177 BeginFastIndex2(WorkN, Work()(Time+TStep,Level,c).bbox(),
00178 Work()(Time+TStep,Level,c).data(), DataType);
00179 DataType uy, vx;
00180 DCoords dx = base::GH().worldStep(ss);
00181 for_2 (i, j, bbi, ss)
00182 vx = (FastIndex2(U,i+ss(0),j)(iv)/FastIndex2(WorkN,i+ss(0),j)-
00183 FastIndex2(U,i-ss(0),j)(iv)/FastIndex2(WorkN,i-ss(0),j))/(2.0*dx(0));
00184 uy = (FastIndex2(U,i,j+ss(1))(iu)/FastIndex2(WorkN,i,j+ss(1))-
00185 FastIndex2(U,i,j-ss(1))(iu)/FastIndex2(WorkN,i,j-ss(1)))/(2.0*dx(1));
00186 FastIndex2(Work,i,j) = std::fabs(vx-uy);
00187 end_for
00188 EndFastIndex2(U);
00189 EndFastIndex2(Work);
00190 EndFastIndex2(WorkN);
00191 end_forall
00192 }
00193 else {
00194 int ipress = base::Dim()+4;
00195 ((output_type*) base::_FileOutput)->Transform(base::U(), base::Work(), Time, Level,
00196 ipress, base::t[Level]);
00197 }
00198
00199 _Accumulation->Accumulate(base::Work(), Time, Level, base::t[Level]);
00200 }
00201
00202 virtual double Tick(int VariableTimeStepping, const double dtv[],
00203 const double cflv[], int& Rejections) {
00204
00205 double cfl = base::Tick(VariableTimeStepping,dtv,cflv,Rejections);
00206 if (!TrackFront) return cfl;
00207
00208 int me = MY_PROC;
00209 int Npoints, finelevel = FineLevel(base::GH());
00210 Npoints = base::shape[0]*RefinedBy(base::GH(),finelevel);
00211 double dx = DeltaX(base::GH(),0,finelevel);
00212
00213 DataType* p = new DataType[Npoints];
00214 DataType* r = new DataType[Npoints];
00215 DataType* xr = new DataType[Npoints];
00216 point_type* xc = new point_type[Npoints];
00217 register int n;
00218
00219 for (n=0; n<Npoints; n++) {
00220 r[n] = base::geom[0]+(n+0.5)*dx;
00221 xc[n](0) = r[n]; xc[n](1) = TrackYPos;
00222 }
00223
00224 for (register int l=0; l<=FineLevel(base::GH()); l++) {
00225 int Time = CurrentTime(base::GH(),l);
00226 int press_cnt = base::Dim()+4;
00227 ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time, l,
00228 press_cnt, base::t[l]);
00229 }
00230
00231 _Interpolation->PointsValues(base::Work(),base::t[0],Npoints,xc,p);
00232 _Interpolation->ArrayCombine(VizServer,Npoints,p);
00233
00234 if (me == VizServer) {
00235 int nc;
00236 f_track(Npoints,p,r,xr,nc,TrackMin,TrackMax,TrackdfMin);
00237 std::ofstream outfile;
00238 std::ostream* out;
00239 std::string fname = "xpos.dat";
00240 outfile.open(fname.c_str(), std::ios::out | std::ios::app);
00241 out = new std::ostream(outfile.rdbuf());
00242 *out << base::t[0] << " " << xr[nc-1] << std::endl;
00243 outfile.close();
00244 delete out;
00245 }
00246
00247 delete [] p;
00248 delete [] r;
00249 delete [] xr;
00250 delete [] xc;
00251
00252 return cfl;
00253 }
00254
00255 protected:
00256 int MaxType, TrackFront;
00257 double TrackYPos, TrackMin, TrackMax, TrackdfMin;
00258 accumulation_type* _Accumulation;
00259 interpolation_type* _Interpolation;