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