00001
00002
00003
00004
00005
00006 #ifndef AMROC_PROBLEM_H
00007 #define AMROC_PROBLEM_H
00008
00009 #include "eulerrhok2.h"
00010
00011 #define NEQUATIONS 14
00012 #define NEQUSED 13
00013 #define NFIXUP 12
00014 #define NAUX 4
00015
00016 #include "ClpProblem.h"
00017
00018 #define OWN_FLAGGING
00019 #define OWN_AMRSOLVER
00020 #include "ClpStdProblem.h"
00021 #include "F77Interfaces/F77UpdateLevelTransfer.h"
00022 #include "SpecialOutput/AMRGridAccumulation.h"
00023
00024 class FlaggingSpecific :
00025 public AMRFlagging<VectorType,FixupType,FlagType,DIM> {
00026 typedef AMRFlagging<VectorType,FixupType,FlagType,DIM> base;
00027 public:
00028 FlaggingSpecific(base::solver_type& solver) : base(solver) {
00029 base::AddCriterion(new ScaledGradient<VectorType,FlagType,DIM>());
00030 base::AddCriterion(new LimiterType<VectorType,FlagType,DIM>());
00031 base::AddCriterion(new AbsoluteError<VectorType,FixupType,FlagType,DIM>(solver));
00032 base::AddCriterion(new RelativeError<VectorType,FixupType,FlagType,DIM>(solver));
00033 base::AddCriterion(new F77ScaledGradient<VectorType,FlagType,DIM>(f_flg));
00034 base::AddCriterion(new F77LimiterType<VectorType,FlagType,DIM>(f_flg));
00035 base::AddCriterion(new F77AbsoluteError<VectorType,FixupType,FlagType,DIM>(solver,f_flg));
00036 base::AddCriterion(new F77RelativeError<VectorType,FixupType,FlagType,DIM>(solver,f_flg));
00037 }
00038 ~FlaggingSpecific() { DeleteAllCriterions(); }
00039 };
00040
00041 class SolverSpecific :
00042 public AMRSolver<VectorType,FixupType,FlagType,DIM> {
00043 typedef VectorType::InternalDataType DataType;
00044 typedef AMRSolver<VectorType,FixupType,FlagType,DIM> base;
00045 typedef AMRGridAccumulation<DataType,DIM> accumulation_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) : base(integ, init, bc) {
00051 SetLevelTransfer(new F77UpdateLevelTransfer<VectorType,DIM>(f_prolong, f_restrict, f_tupdate));
00052 SetFileOutput(new F77FileOutput<VectorType,DIM>(f_out));
00053 SetFixup(new FixupSpecific(integ));
00054 SetFlagging(new FlaggingSpecific(*this));
00055 _Accumulation = new accumulation_type();
00056 MaxType = 1;
00057 }
00058
00059 ~SolverSpecific() {
00060 delete _LevelTransfer;
00061 delete _Flagging;
00062 delete _Fixup;
00063 delete _FileOutput;
00064 delete _Accumulation;
00065 }
00066
00067 virtual void register_at(ControlDevice& Ctrl) { base::register_at(Ctrl); }
00068 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00069 base::register_at(Ctrl, prefix);
00070 RegisterAt(base::LocCtrl,"MaxType",MaxType);
00071 _Accumulation->register_at(base::LocCtrl, prefix);
00072 }
00073
00074 virtual void SetupData() {
00075 base::SetupData();
00076 _Accumulation->SetupData(base::PGH(), base::NGhosts());
00077 }
00078
00079 virtual void Initialize_(const double& dt_start) {
00080 base::Initialize_(dt_start);
00081 _Accumulation->Initialize();
00082 }
00083
00084 virtual void Output() {
00085 int Time = CurrentTime(base::GH(),0);
00086 if (Time == base::LastOutputTime()) return;
00087 _Accumulation->GridCombine(VizServer);
00088 int me = MY_PROC;
00089 if (me == VizServer) {
00090 int Level = _Accumulation->AccLevel();
00091 base::FileOutput_().WriteOut(_Accumulation->AccGrid(),_Accumulation->AccName(),
00092 Time,Level,base::t[0]);
00093 }
00094 base::Output();
00095 }
00096
00097 virtual void Checkpointing_(const char* CheckpointFile) {
00098 base::Checkpointing_(CheckpointFile);
00099 _Accumulation->Checkpointing();
00100 }
00101
00102 virtual bool Restart_(const char* CheckpointFile) {
00103 if (base::Restart_(CheckpointFile)) {
00104 _Accumulation->Restart();
00105 return true;
00106 }
00107 else
00108 return false;
00109 }
00110
00111 virtual void AdvanceLevel(const int Level, int RegridEvery, bool RegridDone,
00112 bool ShadowAllowed, bool DoFixup,
00113 bool RecomposeBaseLev, bool RecomposeHighLev) {
00114
00115 base::AdvanceLevel(Level,RegridEvery,RegridDone,ShadowAllowed,
00116 DoFixup,RecomposeBaseLev,RecomposeHighLev);
00117
00118 int Time = CurrentTime(base::GH(),Level);
00119 if (MaxType) {
00120 AdaptiveBoundaryUpdate(base::U(),Time,Level);
00121 Sync(base::U(),Time,Level);
00122 ExternalBoundaryUpdate(base::U(),Time,Level);
00123
00124 int irho = 1;
00125 int iu = NEQUSED-DIM-2;
00126 int iv = NEQUSED-DIM-1;
00127 int TStep = TimeStep(base::U(),Level);
00128 ((output_type*) base::_FileOutput)->Transform(base::U(), base::Work(), Time+TStep,
00129 Level, irho, base::t[Level]);
00130 forall (base::U(),Time,Level,c)
00131 Coords ss(base::U()(Time,Level,c).bbox().stepsize());
00132 BBox bbi(base::U().interiorbbox(Time,Level,c));
00133 BeginFastIndex2(U, base::U()(Time,Level,c).bbox(),
00134 base::U()(Time,Level,c).data(), VectorType);
00135 BeginFastIndex2(Work, Work()(Time,Level,c).bbox(),
00136 Work()(Time,Level,c).data(), DataType);
00137 BeginFastIndex2(WorkN, Work()(Time+TStep,Level,c).bbox(),
00138 Work()(Time+TStep,Level,c).data(), DataType);
00139 DataType uy, vx;
00140 DCoords dx = base::GH().worldStep(ss);
00141 for_2 (i, j, bbi, ss)
00142 vx = (FastIndex2(U,i+ss(0),j)(iv)/FastIndex2(WorkN,i+ss(0),j)-
00143 FastIndex2(U,i-ss(0),j)(iv)/FastIndex2(WorkN,i-ss(0),j))/(2.0*dx(0));
00144 uy = (FastIndex2(U,i,j+ss(1))(iu)/FastIndex2(WorkN,i,j+ss(1))-
00145 FastIndex2(U,i,j-ss(1))(iu)/FastIndex2(WorkN,i,j-ss(1)))/(2.0*dx(1));
00146 FastIndex2(Work,i,j) = std::fabs(vx-uy);
00147 end_for
00148 EndFastIndex2(U);
00149 EndFastIndex2(Work);
00150 EndFastIndex2(WorkN);
00151 end_forall
00152 }
00153 else {
00154 int ipress = base::Dim()+4;
00155 ((output_type*) base::_FileOutput)->Transform(base::U(), base::Work(), Time, Level,
00156 ipress, base::t[Level]);
00157 }
00158
00159 _Accumulation->Accumulate(base::Work(), Time, Level, base::t[Level]);
00160 }
00161
00162 protected:
00163 int MaxType;
00164 accumulation_type* _Accumulation;
00165 };
00166
00167 #endif
00168
00169