00001
00002
00003
00004
00005
00006 #ifndef AMROC_PROBLEM_H
00007 #define AMROC_PROBLEM_H
00008
00009 #include "eulerm2.h"
00010 #include "ClpProblem.h"
00011
00012 #define f_track FORTRAN_NAME(track, TRACK)
00013 #define f_lsetrfl FORTRAN_NAME(lsrfl, LSRFL)
00014 #define f_lsetrdi FORTRAN_NAME(lsrdi, LSRDI)
00015
00016 extern "C" {
00017 void f_track(const INTEGER &, const DOUBLE *, const DOUBLE *, DOUBLE *, INTEGER&,
00018 const DOUBLE&, const DOUBLE&, const DOUBLE&);
00019 void f_lsetrfl();
00020 void f_lsetrdi();
00021 }
00022
00023 #define OWN_GFMAMRSOLVER
00024 #include "ClpStdGFMProblem.h"
00025 #include "AMRGFMInterpolation.h"
00026
00027 class SolverSpecific :
00028 public AMRGFMSolver<VectorType,FixupType,FlagType,DIM> {
00029 typedef VectorType::InternalDataType DataType;
00030 typedef AMRGFMSolver<VectorType,FixupType,FlagType,DIM> base;
00031 typedef AMRGFMInterpolation<VectorType,FixupType,FlagType,DIM> interpolation_type;
00032 typedef interpolation_type::point_type point_type;
00033 typedef F77GFMFileOutput<VectorType,FixupType,FlagType,DIM> output_type;
00034 public:
00035 SolverSpecific(IntegratorSpecific& integ,
00036 base::initial_condition_type& init,
00037 base::boundary_conditions_type& bc) : base(integ, init, bc) {
00038 SetLevelTransfer(new F77LevelTransfer<VectorType,DIM>(f_prolong, f_restrict));
00039 #ifdef f_flgout
00040 SetFileOutput(new F77GFMFileOutput<VectorType,FixupType,FlagType,DIM>(*this,f_flgout));
00041 #else
00042 SetFileOutput(new GFMFileOutput<VectorType,FixupType,FlagType,DIM>(*this));
00043 #endif
00044 SetFixup(new FixupSpecific(integ));
00045 SetFlagging(new FlaggingSpecific(*this));
00046 AddGFM(new GhostFluidMethod<VectorType,DIM>(
00047 new F77GFMBoundary<VectorType,DIM>(f_ibndex,f_itrans),
00048 new F77GFMLevelSet<DataType,DIM>(f_lset)));
00049 AddGFM(new GhostFluidMethod<VectorType,DIM>(
00050 new F77GFMBoundary<VectorType,DIM>(f_ibndrfl,f_itrans),
00051 new F77GFMLevelSet<DataType,DIM>(f_lsetrdi)));
00052 AddGFM(new GhostFluidMethod<VectorType,DIM>(
00053 new F77GFMBoundary<VectorType,DIM>(f_ibndrfl,f_itrans),
00054 new F77GFMLevelSet<DataType,DIM>(f_lsetrfl)));
00055 _Interpolation = new interpolation_type(*this);
00056 }
00057
00058 ~SolverSpecific() {
00059 DeleteGFM(_GFM[0]);
00060 DeleteGFM(_GFM[1]);
00061 DeleteGFM(_GFM[2]);
00062 delete _Flagging;
00063 delete _Fixup;
00064 delete _Interpolation;
00065 }
00066
00067 virtual void SetupData() {
00068 base::SetupData();
00069 _Interpolation->SetupData(base::PGH(), base::NGhosts());
00070 }
00071
00072 virtual double Tick(int VariableTimeStepping, const double dtv[], const double cflv[],
00073 int& Rejections) {
00074
00075 double cfl = base::Tick(VariableTimeStepping, dtv, cflv, Rejections);
00076
00077 int me = MY_PROC;
00078 int Npoints = 100, npos;
00079 double pi = 4.0*std::atan(1.0), alpha = 0.0, pmax = 0.0;
00080 double shk_pos[10];
00081 double fmin = 0.1, fmax = 1.e4, dfmin = 10.;
00082 for (npos=0; npos<10; alpha+=5.0, npos++) {
00083
00084 DataType* p = new DataType[Npoints];
00085 DataType* r = new DataType[Npoints];
00086 DataType* xr = new DataType[Npoints];
00087 point_type* xc = new point_type[Npoints];
00088 register int n;
00089
00090 for (n=0; n<Npoints; n++) {
00091 r[n] = 8.0/Npoints*n;
00092 xc[n](0) = r[n]*std::cos(alpha*pi/180.); xc[n](1) = r[n]*std::sin(alpha*pi/180.);
00093 }
00094
00095 for (register int l=0; l<=FineLevel(base::GH()); l++) {
00096 int Time = CurrentTime(base::GH(),l);
00097 int press_cnt = base::Dim()+4;
00098 ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time, l,
00099 press_cnt, base::t[l]);
00100 }
00101
00102 _Interpolation->PointsValues(base::Work(),base::t[0],Npoints,xc,p);
00103 _Interpolation->ArrayCombine(VizServer,Npoints,p);
00104
00105 if (me == VizServer) {
00106 int nc;
00107 for (n=0; n<Npoints; n++)
00108 if (p[n]>pmax) pmax = p[n];
00109 if (base::t[0]<0.28) {
00110 dfmin = 10.;
00111 f_track(Npoints,p,r,xr,nc,fmin,fmax,dfmin);
00112 shk_pos[npos] = xr[0];
00113 }
00114 else {
00115 dfmin = 1000.;
00116 f_track(Npoints,p,r,xr,nc,fmin,fmax,dfmin);
00117 shk_pos[npos] = xr[nc-1];
00118 }
00119 }
00120
00121 delete [] p;
00122 delete [] r;
00123 delete [] xr;
00124 delete [] xc;
00125 }
00126
00127 if (me == VizServer) {
00128 std::ofstream outfile;
00129 std::ostream* out;
00130 std::string fname = "shk.dat";
00131 outfile.open(fname.c_str(), std::ios::out | std::ios::app);
00132 out = new std::ostream(outfile.rdbuf());
00133 double sum = 0.0;
00134 for (npos=0; npos<10; npos++)
00135 sum += shk_pos[npos];
00136 *out << base::t[0] << " " << pmax << " " << sum/10.;
00137 for (npos=0; npos<10; npos++)
00138 *out << " " << shk_pos[npos];
00139 *out << std::endl;
00140 outfile.close();
00141 delete out;
00142 }
00143
00144 return cfl;
00145 }
00146
00147 protected:
00148 interpolation_type* _Interpolation;