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