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) {}
00023
00024 virtual void SetGrid(vec_grid_data_type& vec, grid_data_type& workvec, const int& level) {
00025 int Nx = vec.extents(0), Ny = vec.extents(1);
00026 int mxs = base::NGhosts(), mys = base::NGhosts();
00027 int mxe = Nx-base::NGhosts()-1, mye = Ny-base::NGhosts()-1;
00028 VectorType *mesh = (VectorType *)vec.databuffer();
00029
00030 DCoords lbcorner = base::GH().worldCoords(vec.lower(), vec.stepsize());
00031 DCoords dx = base::GH().worldStep(vec.stepsize());
00032
00033 for (register int j=mys; j<=mye; j++) {
00034 double yc = (j+0.5)*dx(1)+lbcorner(1);
00035 for (register int i=mxs; i<=mxe; i++) {
00036 double xc = (i+0.5)*dx(0)+lbcorner(0);
00037 register int pos = Scheme().idx(i,j,Nx);
00038
00039 mesh[pos](0) = Scheme().Gamma()*Scheme().Gamma();
00040 mesh[pos](1) = Scheme().Gamma();
00041 mesh[pos](2) = 0.;
00042 mesh[pos](3) = -std::sin(yc);
00043 mesh[pos](4) = std::sin(xc);
00044 mesh[pos](5) = 0.;
00045 mesh[pos](6) = -std::sin(yc);
00046 mesh[pos](7) = std::sin(2.*xc);
00047 mesh[pos](8) = 0.;
00048 }
00049 }
00050
00051 grid_data_type rhoE(vec.bbox());
00052 Scheme().computeRhoE(vec,rhoE);
00053 Scheme().primitive2conservative(vec,rhoE);
00054 }
00055 };
00056
00057 class SolverSpecific :
00058 public AMRSolver<VectorType,FixupType,FlagType,DIM> {
00059 typedef AMRSolver<VectorType,FixupType,FlagType,DIM> base;
00060 typedef VectorType::InternalDataType DataType;
00061 typedef SchemeFileOutput<SchemeType,DIM> output_type;
00062 public:
00063 SolverSpecific(IntegratorSpecific& integ,
00064 base::initial_condition_type& init,
00065 base::boundary_conditions_type& bc) : base(integ, init, bc),
00066 DivergenceOutput(1), ConservationOutput(1) {
00067 SetLevelTransfer(new F77LevelTransfer<VectorType,DIM>(f_prolong, f_restrict));
00068 SetFileOutput(new SchemeFileOutput<SchemeType,DIM>(integ.Scheme()));
00069 SetFixup(new FixupSpecific());
00070 SetFlagging(new FlaggingSpecific(*this));
00071 }
00072
00073 ~SolverSpecific() {
00074 delete _LevelTransfer;
00075 delete _Flagging;
00076 delete _Fixup;
00077 delete _FileOutput;
00078 }
00079
00080 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00081 base::register_at(Ctrl,prefix);
00082 RegisterAt(base::LocCtrl,"DivergenceOutput",DivergenceOutput);
00083 RegisterAt(base::LocCtrl,"ConservationOutput",ConservationOutput);
00084 }
00085
00086 virtual void AfterLevelStep(const int Level) {
00087 if (DivergenceOutput==1) {
00088 int Time = CurrentTime(base::GH(),Level);
00089 int idiv = 12;
00090 Work().GF_Fill(DataType(0.),Time,Level);
00091 ((output_type*) base::_FileOutput)->Transform(base::U(), base::Work(), Time,
00092 Level, idiv, base::t[Level]);
00093 DCoords dx = base::GH().worldStep(Work().GF_StepSize(Level));
00094 DataType divSum = Work().GF_sum(Time,Level)*(dx(0)*dx(1));
00095 DataType divMax = Work().GF_norm(Time,Level,DAGHNormMax);
00096 DataType divSumAbs = Work().GF_norm(Time,Level,DAGHNormL1);
00097 if (MY_PROC==VizServer) {
00098 std::ofstream outfile;
00099 std::ostream* out;
00100 std::string fname = "divergence.txt";
00101 outfile.open(fname.c_str(), std::ios::out | std::ios::app);
00102 out = new std::ostream(outfile.rdbuf());
00103 *out << base::t[Level] << ": \t After Level=" << Level << ": \t|divMax|= " << divMax << " \t Int div(i,j)= "
00104 << divSum <<" \t Int |div(i,j)|= " << divSumAbs << std::endl;
00105 outfile.close();
00106 }
00107 }
00108 }
00109
00110 virtual double Tick(int VariableTimeStepping, const double dtv[],
00111 const double cflv[], int& Rejections) {
00112 double cfl = base::Tick(VariableTimeStepping, dtv, cflv, Rejections);
00113 if (ConservationOutput==1) {
00114 DataType E = IntegrateOutputQuantity(4,0);
00115 DataType div = IntegrateOutputQuantity(12,0);
00116 DataType dmax = IntegrateOutputQuantity(12,1);
00117 DataType H = IntegrateOutputQuantity(16,0);
00118 if (MY_PROC==VizServer) {
00119 std::ofstream outfile;
00120 std::ostream* out;
00121 std::string fname = "conservation.txt";
00122 outfile.open(fname.c_str(), std::ios::out | std::ios::app);
00123 out = new std::ostream(outfile.rdbuf());
00124 *out << base::t[0] << " " << E << " " << div << " " << dmax << " " << H << std::endl;
00125 outfile.close();
00126 }
00127 }
00128 return cfl;
00129 }
00130
00131 double IntegrateOutputQuantity(int idiv, int param) {
00132 DataType ret = 0.;
00133 int Time = CurrentTime(base::GH(),0);
00134 for (int lev=0; lev<=FineLevel(base::GH()); lev++) {
00135 DCoords dx = base::GH().worldStep(base::Work().GF_StepSize(lev));
00136 base::Work().GF_Fill(DataType(0.),Time,lev);
00137 ((output_type*) base::_FileOutput)->Transform(base::U(), base::Work(), Time,
00138 lev, idiv, base::t[lev]);
00139 forall (base::Work(),Time,lev,c)
00140 if (lev<FineLevel(base::GH())) {
00141 forall (base::Work(),Time,lev+1,c2)
00142 BBox bb(coarsen(base::Work().interiorbbox(Time,lev+1,c2),base::GH().refinefactor(lev)));
00143 base::Work()(Time,lev,c).equals(DataType(0.0),bb);
00144 end_forall
00145 }
00146 end_forall
00147 if (param==0)
00148 ret += base::Work().GF_sum(Time,lev)*(dx(0)*dx(1));
00149 else
00150 ret = std::max(ret,base::Work().GF_norm(Time,lev,DAGHNormMax));
00151 }
00152 return ret;
00153 }
00154
00155 private:
00156 int DivergenceOutput, ConservationOutput;
00157 };