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