00001
00002
00003
00004
00005
00006 #ifndef AMROC_PROBLEM_H
00007 #define AMROC_PROBLEM_H
00008
00009 #include "eulerrhok2.h"
00010 #define NEQUATIONS 6
00011 #define NEQUSED 6
00012 #define NFIXUP 5
00013 #define NAUX 1
00014
00015 #include "ClpProblem.h"
00016
00017 #define OWN_FLAGGING
00018 #define OWN_AMRSOLVER
00019 #include "ClpStdProblem.h"
00020 #include "F77Interfaces/F77UpdateLevelTransfer.h"
00021 #include "AMRInterpolation.h"
00022
00023 class FlaggingSpecific :
00024 public AMRFlagging<VectorType,FixupType,FlagType,DIM> {
00025 typedef AMRFlagging<VectorType,FixupType,FlagType,DIM> base;
00026 public:
00027 FlaggingSpecific(base::solver_type& solver) : base(solver) {
00028 base::AddCriterion(new ScaledGradient<VectorType,FlagType,DIM>());
00029 base::AddCriterion(new LimiterType<VectorType,FlagType,DIM>());
00030 base::AddCriterion(new AbsoluteError<VectorType,FixupType,FlagType,DIM>(solver));
00031 base::AddCriterion(new RelativeError<VectorType,FixupType,FlagType,DIM>(solver));
00032 base::AddCriterion(new F77ScaledGradient<VectorType,FlagType,DIM>(f_flg));
00033 base::AddCriterion(new F77LimiterType<VectorType,FlagType,DIM>(f_flg));
00034 base::AddCriterion(new F77AbsoluteError<VectorType,FixupType,FlagType,DIM>(solver,f_flg));
00035 base::AddCriterion(new F77RelativeError<VectorType,FixupType,FlagType,DIM>(solver,f_flg));
00036 }
00037 ~FlaggingSpecific() { DeleteAllCriterions(); }
00038 };
00039
00040 class SolverSpecific :
00041 public AMRSolver<VectorType,FixupType,FlagType,DIM> {
00042 typedef AMRSolver<VectorType,FixupType,FlagType,DIM> base;
00043 typedef VectorType::InternalDataType DataType;
00044 typedef AMRInterpolation<VectorType,DIM> interpolation_type;
00045 typedef interpolation_type::point_type point_type;
00046 typedef F77FileOutput<VectorType,DIM> output_type;
00047 public:
00048 SolverSpecific(IntegratorSpecific& integ,
00049 base::initial_condition_type& init,
00050 base::boundary_conditions_type& bc) :
00051 AMRSolver<VectorType,FixupType,FlagType,DIM>(integ, init, bc), OutputPressureEvery(1) {
00052 SetLevelTransfer(new F77UpdateLevelTransfer<VectorType,DIM>(f_prolong, f_restrict, f_tupdate));
00053 SetFileOutput(new F77FileOutput<VectorType,DIM>(f_out));
00054 SetFixup(new FixupSpecific(integ));
00055 SetFlagging(new FlaggingSpecific(*this));
00056 _Interpolation = new interpolation_type();
00057 }
00058
00059 ~SolverSpecific() {
00060 delete _LevelTransfer;
00061 delete _Flagging;
00062 delete _Fixup;
00063 delete _FileOutput;
00064 delete _Interpolation;
00065 }
00066
00067 virtual void SetupData() {
00068 base::SetupData();
00069 _Interpolation->SetupData(base::PGH(), base::NGhosts());
00070 }
00071 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00072 base::register_at(Ctrl,prefix);
00073 RegisterAt(base::LocCtrl,"OutputPressureEvery",OutputPressureEvery);
00074 }
00075 virtual void register_at(ControlDevice& Ctrl) {
00076 base::register_at(Ctrl);
00077 }
00078
00079
00080
00081 virtual double Tick(int VariableTimeStepping, const double dtv[],
00082 const double cflv[], int& Rejections) {
00083
00084 double cfl = base::Tick(VariableTimeStepping, dtv, cflv, Rejections);
00085
00086 int TimeZero = CurrentTime(base::GH(),0);
00087 TimeZero /= RefinedBy(base::GH(),MaxLevel(base::GH()));
00088 if (TimeZero%OutputPressureEvery!=0 ||
00089 base::GH().nbndry()<=1) return cfl;
00090
00091 int me = MY_PROC;
00092
00093
00094 for (register int l=0; l<=FineLevel(base::GH()); l++) {
00095 int Time = CurrentTime(base::GH(),l);
00096 int press_cnt = base::Dim()+4;
00097 ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time, l,
00098 press_cnt, base::t[l]);
00099 }
00100
00101 const double* bndrycoords = base::GH().wholebndry();
00102 double lpy = bndrycoords[1+base::GH().nbndry()*(0+2*1)];
00103 double upy = bndrycoords[1+base::GH().nbndry()*(1+2*1)];
00104 double upx = bndrycoords[1+base::GH().nbndry()*(1+2*0)];
00105 double dy = base::GH().delta_x(1,MaxLevel(base::GH()));
00106 double pav = 0., global_pav = 0.;
00107 int Npoints = int((upy-lpy-0.5*dy)/dy);
00108
00109 DataType* p = new DataType[Npoints];
00110 point_type* xc = new point_type[Npoints];
00111 register int n;
00112 for (n=0; n<Npoints; n++) {
00113 xc[n](0) = upx;
00114 xc[n](1) = lpy+(n+0.5)*dy;
00115 }
00116
00117 _Interpolation->PointsValues(base::Work(),base::t[0],Npoints,xc,p);
00118 for (n=0; n<Npoints; n++)
00119 if (p[n]>-1.e37)
00120 pav += p[n];
00121 delete [] p;
00122 delete [] xc;
00123
00124 #ifdef DAGH_NO_MPI
00125 global_pav = pav;
00126 #else
00127 int Nval = 1;
00128 MPI_Allreduce(&pav, &global_pav, Nval, MPI_DOUBLE, MPI_SUM, base::Comm());
00129 #endif
00130
00131
00132 if (me == VizServer) {
00133 std::ofstream outfile;
00134 std::ostream* out;
00135 std::string fname = "pressure.dat";
00136 outfile.open(fname.c_str(), std::ios::out | std::ios::app);
00137 out = new std::ostream(outfile.rdbuf());
00138 *out << base::t[0] << " " << global_pav/Npoints << std::endl;
00139 outfile.close();
00140 delete out;
00141 }
00142
00143 return cfl;
00144 }
00145
00146 protected:
00147 int OutputPressureEvery;
00148 interpolation_type* _Interpolation;
00149 };
00150