00001
00002
00003 #ifndef AMROC_PROBLEM_H
00004 #define AMROC_PROBLEM_H
00005
00006 #define DIM 1
00007
00008
00009 #define OWN_AMRSOLVER
00010 #include "WENOProblem.h"
00011
00012 #define f_exact FORTRAN_NAME(exact, EXACT)
00013
00014 extern "C" {
00015 void f_exact();
00016 }
00017
00018 #ifdef SYNCTIME
00019 #include "WENOStdSyncTimeProblem.h"
00020 #else
00021 #include "WENOStdProblem.h"
00022 #endif
00023 #include "F77Interfaces/F77ExactSolution.h"
00024
00025 class SolverSpecific :
00026 #ifdef SYNCTIME
00027 public AMRPreAdaptSolverSyncTime<VectorType,FixupType,FlagType,DIM> {
00028 typedef AMRPreAdaptSolverSyncTime<VectorType,FixupType,FlagType,DIM> base;
00029 #else
00030 public AMRPreAdaptSolver<VectorType,FixupType,FlagType,DIM> {
00031 typedef AMRPreAdaptSolver<VectorType,FixupType,FlagType,DIM> base;
00032 #endif
00033 typedef AMRSolver<VectorType,FixupType,FlagType,DIM> bbase;
00034 typedef WENOIntegrator<VectorType,DIM> weno_integ_type;
00035 typedef VectorType::InternalDataType DataType;
00036 public:
00037 SolverSpecific(IntegratorSpecific& integ,
00038 bbase::initial_condition_type& init,
00039 bbase::boundary_conditions_type& bc) :
00040 #ifdef SYNCTIME
00041 AMRPreAdaptSolverSyncTime<VectorType,FixupType,FlagType,DIM>(integ, init, bc) {
00042 #else
00043 AMRPreAdaptSolver<VectorType,FixupType,FlagType,DIM>(integ, init, bc) {
00044 #endif
00045 SetLevelTransfer(new F77LevelTransfer<VectorType,DIM>(f_prolong, f_restrict));
00046 #ifdef f_flgout
00047 SetFileOutput(new F77FileOutput<VectorType,DIM>(f_flgout));
00048 #else
00049 SetFileOutput(new FileOutput<VectorType,DIM>());
00050 #endif
00051 SetFixup(new FixupSpecific());
00052 SetFlagging(new FlaggingSpecific(*this));
00053
00054 SetExactSolution(new F77ExactSolution<VectorType,DIM>(f_exact));
00055 }
00056
00057 ~SolverSpecific() {
00058 delete _LevelTransfer;
00059 delete _Flagging;
00060 delete _Fixup;
00061 delete _FileOutput;
00062 }
00063
00064 virtual void Initialize_(const double& dt_start) {
00065 base::Initialize_(dt_start);
00066 CalculateCurrentMass(mass_old);
00067 }
00068
00069 virtual double Tick(int VariableTimeStepping, const double dtv[],
00070 const double cflv[], int& Rejections) {
00071
00072 double cfl = base::Tick(VariableTimeStepping,dtv,cflv,Rejections);
00073
00074 VectorType mass_new;
00075 CalculateCurrentMass(mass_new);
00076
00077
00078 int me = MY_PROC;
00079 if (me == VizServer) {
00080 std::ofstream outfile;
00081 std::ostream* out;
00082 std::string fname = "mass.dat";
00083 outfile.open(fname.c_str(), std::ios::out | std::ios::app);
00084 out = new std::ostream(outfile.rdbuf());
00085 *out << base::t[0];
00086 for (int n=0; n<base::NEquations(); n++)
00087 *out << " " << mass_old[n]-mass_new[n];
00088 *out << std::endl;
00089 outfile.close();
00090 delete out;
00091 }
00092
00093
00094 int irho = 0;
00095 int iu = 1;
00096 int iv = 2;
00097 double kin = 0.0;
00098 for ( int Level = FineLevel(base::GH()) ; Level >= 0 ; Level --)
00099 {
00100 int Time = CurrentTime(base::GH(),Level);
00101 forall (base::U(),Time,Level,c)
00102 Coords ss(base::U()(Time,Level,c).bbox().stepsize());
00103 BBox bbi(base::U().interiorbbox(Time,Level,c));
00104 BeginFastIndex2(U, base::U()(Time,Level,c).bbox(),
00105 base::U()(Time,Level,c).data(), VectorType);
00106 BeginFastIndex2(Work, Work()(Time,Level,c).bbox(),
00107 Work()(Time,Level,c).data(), DataType);
00108 DataType rho, u, v;
00109 for_2 (i, j, bbi, ss)
00110 rho = FastIndex2(U,i,j)(irho);
00111 u = FastIndex2(U,i,j)(iu)/rho;
00112 v = FastIndex2(U,i,j)(iv)/rho;
00113 FastIndex2(Work,i,j) = (u*u+v*v)/2.0;
00114 end_for
00115 EndFastIndex2(U);
00116 EndFastIndex2(Work);
00117 end_forall
00118
00119 forall (base::Work(),Time,Level,c)
00120 if (Level<FineLevel(base::GH())) {
00121 forall (base::Work(),Time,Level+1,c2)
00122 BBox bb(coarsen(base::Work().interiorbbox(Time,Level+1,c2),base::GH().refinefactor(Level)));
00123 base::Work()(Time,Level,c).equals(DataType(0.0),bb);
00124 end_forall
00125 }
00126 end_forall
00127 kin += base::Work().GF_sum(Time,Level)/(Level+1)/(Level+1);
00128 }
00129 DCoords dx = base::GH().worldStep(base::U().GF_StepSize(0));
00130 for (int d=0; d<base::Dim(); d++) kin *= dx(d);
00131
00132 if (MY_PROC == VizServer) {
00133 std::ofstream outfile;
00134 std::ostream* out;
00135 std::string fname = "stats.dat";
00136 outfile.open(fname.c_str(), std::ios::out | std::ios::app);
00137 out = new std::ostream(outfile.rdbuf());
00138 out->setf(std::ios::scientific, std::ios::floatfield);
00139 *out << base::t[0] << " " << kin << std::endl;
00140 outfile.close();
00141 delete out;
00142 }
00143
00144 return cfl;
00145 }
00146
00147 void CalculateCurrentMass(VectorType& mass) {
00148 int Level=0;
00149 int Time = CurrentTime(base::GH(),Level);
00150 int TStep = TimeStep(base::U(),Level);
00151 forall (base::U(),Time,Level,c)
00152 base::U()(Time+TStep,Level,c).equals(base::U()(Time,Level,c));
00153 end_forall
00154
00155 DCoords dx = base::GH().worldStep(base::U().GF_StepSize(Level));
00156 mass = Sum(base::U(),Time+TStep,Level);
00157 for (int d=0; d<base::Dim(); d++) mass *= dx(d);
00158 }
00159
00160 protected:
00161 VectorType mass_old;
00162 };
00163