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