00001
00002
00003 #ifndef AMROC_PROBLEM_H
00004 #define AMROC_PROBLEM_H
00005
00006 #define DIM 2
00007
00008 #include "MHDProblem.h"
00009 #include "EGLM2D.h"
00010 #include <iomanip>
00011
00012 typedef double DataType;
00013 typedef EGLM2D<DataType> SchemeType;
00014
00015 #define OWN_INITIALCONDITION
00016 #define OWN_AMRSOLVER
00017 #include "MHDStdProblem.h"
00018
00019 class InitialConditionSpecific : public SchemeInitialCondition<SchemeType,DIM> {
00020 typedef SchemeInitialCondition<SchemeType,DIM> base;
00021 public:
00022 InitialConditionSpecific(SchemeType &scheme) : base(scheme), iccase(0) {}
00023
00024 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00025 base::register_at(Ctrl,prefix);
00026 RegisterAt(base::LocCtrl,"IC_case",iccase);
00027 }
00028
00029 virtual void SetGrid(vec_grid_data_type& fvec, grid_data_type& workvec, const int& level) {
00030 if (iccase==1) {
00031 DataType vec[19] = { 1., 20., 0., 10., 0., 0., 5., 5., 0., 1., 1., 0., -10., 0., 0., 5., 5., 0., 0. };
00032 DataType auxKonst = sqrt(DataType(16.)*atan(DataType(1.)));
00033 vec[6] = vec[6]/auxKonst;
00034 vec[7] = vec[7]/auxKonst;
00035 vec[15] = vec[15]/auxKonst;
00036 vec[16] = vec[16]/auxKonst;
00037 Scheme().ICStandard(fvec,SchemeType::DiscontinuityX,vec,19,SchemeType::Primitive);
00038 }
00039 else if (iccase==2) {
00040 DataType vec[19] = { 1.08, 0.95, 0., 1.2, 0.01, 0.5, 2., 3.6, 2., 1., 1., 0., 0., 0., 0., 2., 4., 2., 0. };
00041 DataType auxKonst = sqrt(DataType(16.)*atan(DataType(1.)));
00042 DataType auxKonstLi = sqrt(DataType(8.)*atan(DataType(1.)));
00043 vec[6] = vec[6]/auxKonst;
00044 vec[7] = vec[7]/auxKonst;
00045 vec[8] = vec[8]/auxKonstLi;
00046 vec[15] = vec[15]/auxKonst;
00047 vec[16] = vec[16]/auxKonst;
00048 vec[17] = vec[17]/auxKonstLi;
00049 Scheme().ICStandard(fvec,SchemeType::DiscontinuityX,vec,19,SchemeType::Primitive);
00050 }
00051 else if (iccase==3) {
00052 DataType vec[19] = { 1.08, 0.95, 0., 0.01, 1.2, 0.5, 3.6, 2., 2., 1., 1., 0., 0., 0., 0., 4., 2., 2., 0. };
00053 DataType auxKonst = sqrt(DataType(16.)*atan(DataType(1.)));
00054 DataType auxKonstLi = sqrt(DataType(8.)*atan(DataType(1.)));
00055 vec[6] = vec[6]/auxKonst;
00056 vec[7] = vec[7]/auxKonst;
00057 vec[8] = vec[8]/auxKonstLi;
00058 vec[15] = vec[15]/auxKonst;
00059 vec[16] = vec[16]/auxKonst;
00060 vec[17] = vec[17]/auxKonstLi;
00061 Scheme().ICStandard(fvec,SchemeType::DiscontinuityY,vec,19,SchemeType::Primitive);
00062 }
00063 }
00064
00065 protected:
00066 int iccase;
00067 };
00068
00069 class SolverSpecific :
00070 public AMRPreAdaptSolver<VectorType,FixupType,FlagType,DIM> {
00071 typedef AMRPreAdaptSolver<VectorType,FixupType,FlagType,DIM> base;
00072 typedef VectorType::InternalDataType DataType;
00073 typedef SchemeFileOutput<SchemeType,DIM> output_type;
00074 public:
00075 SolverSpecific(IntegratorSpecific& integ,
00076 base::initial_condition_type& init,
00077 base::boundary_conditions_type& bc) : base(integ, init, bc),
00078 DivergenceOutput(1), ConservationOutput(1) {
00079 SetLevelTransfer(new F77LevelTransfer<VectorType,DIM>(f_prolong, f_restrict));
00080 SetFileOutput(new SchemeFileOutput<SchemeType,DIM>(integ.Scheme()));
00081 SetFixup(new FixupSpecific());
00082 SetFlagging(new FlaggingSpecific(*this));
00083 }
00084
00085 ~SolverSpecific() {
00086 delete _LevelTransfer;
00087 delete _Flagging;
00088 delete _Fixup;
00089 delete _FileOutput;
00090 }
00091
00092 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00093 base::register_at(Ctrl,prefix);
00094 RegisterAt(base::LocCtrl,"DivergenceOutput",DivergenceOutput);
00095 RegisterAt(base::LocCtrl,"ConservationOutput",ConservationOutput);
00096 }
00097
00098 virtual void AfterLevelStep(const int Level) {
00099 if (DivergenceOutput==1) {
00100 int Time = CurrentTime(base::GH(),Level);
00101 int idiv = 12;
00102 Work().GF_Fill(DataType(0.),Time,Level);
00103 ((output_type*) base::_FileOutput)->Transform(base::U(), base::Work(), Time,
00104 Level, idiv, base::t[Level]);
00105 DCoords dx = base::GH().worldStep(Work().GF_StepSize(Level));
00106 DataType divSum = Work().GF_sum(Time,Level)*(dx(0)*dx(1));
00107 DataType divMax = Work().GF_norm(Time,Level,DAGHNormMax);
00108 DataType divSumAbs = Work().GF_norm(Time,Level,DAGHNormL1);
00109 if (MY_PROC==VizServer) {
00110 std::ofstream outfile;
00111 std::ostream* out;
00112 std::string fname = "divergence.txt";
00113 outfile.open(fname.c_str(), std::ios::out | std::ios::app);
00114 out = new std::ostream(outfile.rdbuf());
00115 *out << base::t[Level] << ": \t After Level=" << Level << ": \t|divMax|= " << divMax << " \t Int div(i,j)= "
00116 << divSum <<" \t Int |div(i,j)|= " << divSumAbs << std::endl;
00117 outfile.close();
00118 }
00119 }
00120 }
00121
00122 virtual double Tick(int VariableTimeStepping, const double dtv[],
00123 const double cflv[], int& Rejections) {
00124 double cfl = base::Tick(VariableTimeStepping, dtv, cflv, Rejections);
00125 if (ConservationOutput==1) {
00126 DataType E = IntegrateOutputQuantity(4,0);
00127 DataType div = IntegrateOutputQuantity(12,0);
00128 DataType dmax = IntegrateOutputQuantity(12,1);
00129 DataType H = IntegrateOutputQuantity(16,0);
00130 if (MY_PROC==VizServer) {
00131 std::ofstream outfile;
00132 std::ostream* out;
00133 std::string fname = "conservation.txt";
00134 outfile.open(fname.c_str(), std::ios::out | std::ios::app);
00135 out = new std::ostream(outfile.rdbuf());
00136 *out << base::t[0] << " " << E << " " << div << " " << dmax << " " << H << std::endl;
00137 outfile.close();
00138 }
00139 }
00140 return cfl;
00141 }
00142
00143 double IntegrateOutputQuantity(int idiv, int param) {
00144 DataType ret = 0.;
00145 int Time = CurrentTime(base::GH(),0);
00146 for (int lev=0; lev<=FineLevel(base::GH()); lev++) {
00147 DCoords dx = base::GH().worldStep(base::Work().GF_StepSize(lev));
00148 base::Work().GF_Fill(DataType(0.),Time,lev);
00149 ((output_type*) base::_FileOutput)->Transform(base::U(), base::Work(), Time,
00150 lev, idiv, base::t[lev]);
00151 forall (base::Work(),Time,lev,c)
00152 if (lev<FineLevel(base::GH())) {
00153 forall (base::Work(),Time,lev+1,c2)
00154 BBox bb(coarsen(base::Work().interiorbbox(Time,lev+1,c2),base::GH().refinefactor(lev)));
00155 base::Work()(Time,lev,c).equals(DataType(0.0),bb);
00156 end_forall
00157 }
00158 end_forall
00159 if (param==0)
00160 ret += base::Work().GF_sum(Time,lev)*(dx(0)*dx(1));
00161 else
00162 ret = std::max(ret,base::Work().GF_norm(Time,lev,DAGHNormMax));
00163 }
00164 return ret;
00165 }
00166
00167 private:
00168 int DivergenceOutput, ConservationOutput;
00169 };