00001
00002
00003 #ifndef AMROC_PROBLEM_H
00004 #define AMROC_PROBLEM_H
00005
00006 #define DIM 2
00007
00008
00009 #include "WENOProblem.h"
00010
00011 #define f_exact FORTRAN_NAME(exact, EXACT)
00012
00013 extern "C" {
00014 void f_exact();
00015 }
00016
00017 #define OWN_AMRSOLVER
00018 #ifdef SYNCTIME
00019 #include "WENOStdSyncTimeProblem.h"
00020 #else
00021 #include "WENOStdProblem.h"
00022 #endif
00023 #include "F77Interfaces/F77ExactSolution.h"
00024 class SolverSpecific :
00025 #ifdef SYNCTIME
00026 public AMRSolverSyncTime<VectorType,FixupType,FlagType,DIM> {
00027 typedef AMRSolverSyncTime<VectorType,FixupType,FlagType,DIM> base;
00028 #else
00029 public AMRSolver<VectorType,FixupType,FlagType,DIM> {
00030 typedef AMRSolver<VectorType,FixupType,FlagType,DIM> base;
00031 #endif
00032 typedef WENOFixup<VectorType,FixupType,DIM> weno_fixup_type;
00033 typedef WENOIntegrator<VectorType,DIM> weno_integ_type;
00034 public:
00035 SolverSpecific(IntegratorSpecific& integ,
00036 base::initial_condition_type& init,
00037 base::boundary_conditions_type& bc) :
00038 #ifdef SYNCTIME
00039 AMRSolverSyncTime<VectorType,FixupType,FlagType,DIM>(integ, init, bc) {
00040 #else
00041 AMRSolver<VectorType,FixupType,FlagType,DIM>(integ, init, bc) {
00042 #endif
00043 SetLevelTransfer(new F77LevelTransfer<VectorType,DIM>(f_prolong, f_restrict));
00044 #ifdef f_flgout
00045 SetFileOutput(new F77FileOutput<VectorType,DIM>(f_flgout));
00046 #else
00047 SetFileOutput(new FileOutput<VectorType,DIM>());
00048 #endif
00049 SetFixup(new FixupSpecific());
00050 SetFlagging(new FlaggingSpecific(*this));
00051
00052 SetExactSolution(new F77ExactSolution<VectorType,DIM>(f_exact));
00053 }
00054
00055 ~SolverSpecific() {
00056 delete _LevelTransfer;
00057 delete _Flagging;
00058 delete _Fixup;
00059 delete _FileOutput;
00060 }
00061
00062 virtual void Initialize_(const double& dt_start) {
00063 base::Initialize_(dt_start);
00064 CalculateCurrentMass(mass_old);
00065 }
00066
00067 virtual double Tick(int VariableTimeStepping, const double dtv[],
00068 const double cflv[], int& Rejections) {
00069
00070 double cfl = base::Tick(VariableTimeStepping,dtv,cflv,Rejections);
00071
00072 VectorType mass_new;
00073 CalculateCurrentMass(mass_new);
00074
00075
00076 int me = MY_PROC;
00077 if (me == VizServer) {
00078 std::ofstream outfile;
00079 std::ostream* out;
00080 std::string fname = "mass.dat";
00081 outfile.open(fname.c_str(), std::ios::out | std::ios::app);
00082 out = new std::ostream(outfile.rdbuf());
00083 *out << base::t[0];
00084 for (int n=0; n<base::NEquations(); n++)
00085 *out << " " << mass_old[n]-mass_new[n];
00086 *out << std::endl;
00087 outfile.close();
00088 delete out;
00089 }
00090
00091 return cfl;
00092 }
00093
00094 void CalculateCurrentMass(VectorType& mass) {
00095 int Level=0;
00096 int Time = CurrentTime(base::GH(),Level);
00097 int TStep = TimeStep(base::U(),Level);
00098 forall (base::U(),Time,Level,c)
00099 base::U()(Time+TStep,Level,c).equals(base::U()(Time,Level,c));
00100 end_forall
00101
00102 DCoords dx = base::GH().worldStep(base::U().GF_StepSize(Level));
00103 mass = Sum(base::U(),Time+TStep,Level);
00104 for (int d=0; d<base::Dim(); d++) mass *= dx(d);
00105 }
00106
00107 protected:
00108 VectorType mass_old;