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 #include "AMRGFMSolver.h"
00019 #include "Interfaces/SchemeGFMBoundary.h"
00020 #include "Interfaces/SchemeGFMFileOutput.h"
00021
00022 DataType rho0=1.;
00023 DataType rd, xs[2], rs, T;
00024 DataType pi=4.0*std::atan(1.0);
00025
00026
00027 class InitialConditionSpecific : public SchemeInitialCondition<SchemeType,DIM> {
00028 typedef SchemeInitialCondition<SchemeType,DIM> base;
00029 public:
00030 InitialConditionSpecific(SchemeType &scheme) : base(scheme), iccase(0) {}
00031
00032 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00033 base::register_at(Ctrl,prefix);
00034 RegisterAt(base::LocCtrl,"IC_case",iccase);
00035 }
00036
00037 virtual void SetGrid(vec_grid_data_type& fvec, grid_data_type& workvec, const int& level) {
00038 if (iccase==1) {
00039 DataType vec[19] = { 1., 20., 0., 10., 0., 0., 5., 5., 0., 1., 1., 0., -10., 0., 0., 5., 5., 0., 0. };
00040 DataType auxKonst = sqrt(DataType(16.)*atan(DataType(1.)));
00041 vec[6] = vec[6]/auxKonst;
00042 vec[7] = vec[7]/auxKonst;
00043 vec[15] = vec[15]/auxKonst;
00044 vec[16] = vec[16]/auxKonst;
00045 Scheme().ICStandard(fvec,SchemeType::DiscontinuityX,vec,19,SchemeType::Primitive);
00046 }
00047 else if (iccase==2) {
00048 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. };
00049 DataType auxKonst = sqrt(DataType(16.)*atan(DataType(1.)));
00050 DataType auxKonstLi = sqrt(DataType(8.)*atan(DataType(1.)));
00051 vec[6] = vec[6]/auxKonst;
00052 vec[7] = vec[7]/auxKonst;
00053 vec[8] = vec[8]/auxKonstLi;
00054 vec[15] = vec[15]/auxKonst;
00055 vec[16] = vec[16]/auxKonst;
00056 vec[17] = vec[17]/auxKonstLi;
00057 Scheme().ICStandard(fvec,SchemeType::DiscontinuityX,vec,19,SchemeType::Primitive);
00058 }
00059 else if (iccase==3) {
00060 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. };
00061 DataType auxKonst = sqrt(DataType(16.)*atan(DataType(1.)));
00062 DataType auxKonstLi = sqrt(DataType(8.)*atan(DataType(1.)));
00063 vec[6] = vec[6]/auxKonst;
00064 vec[7] = vec[7]/auxKonst;
00065 vec[8] = vec[8]/auxKonstLi;
00066 vec[15] = vec[15]/auxKonst;
00067 vec[16] = vec[16]/auxKonst;
00068 vec[17] = vec[17]/auxKonstLi;
00069 Scheme().ICStandard(fvec,SchemeType::DiscontinuityY,vec,19,SchemeType::Primitive);
00070 }
00071 }
00072
00073 protected:
00074 int iccase;
00075 };
00076
00077
00078 class GFMBoundarySpecific : public SchemeGFMBoundary<SchemeType,DIM> {
00079 typedef SchemeGFMBoundary<SchemeType,DIM> base;
00080 public:
00081 GFMBoundarySpecific(SchemeType &scheme) : base(scheme) {}
00082
00083 virtual void SetBndry(vec_grid_data_type& gdu, const int& nc, const int* idx,
00084 const VectorType* u, const point_type* xc, const DataType* distance,
00085 const point_type* normal, DataType* aux, const int naux, const double t)
00086 { Scheme().GFMBCStandard(gdu,base::Type,nc,idx,u,xc,distance,normal,aux,naux); }
00087
00088 virtual void SetBndryAux(vec_grid_data_type& gdu, grid_data_type& gdphi,
00089 const VectorType* u, DataType* aux, const int& Level,
00090 double t, const int& nc, const int* idx,
00091 const point_type* xc, const DataType* distance,
00092 const point_type* normal) {
00093 DataType xsc[2], vc[2];
00094 xsc[0] = xs[0]; xsc[1] = xs[1];
00095 vc[0] = 0.; vc[1] = 0.;
00096 if (T!=0.) {
00097 xsc[0] += rs*std::cos(2.*pi/T*t);
00098 xsc[1] += rs*std::sin(2.*pi/T*t);
00099 vc[0] = -2.*pi/T*std::sin(2.*pi/T*t);
00100 vc[1] = 2.*pi/T*std::cos(2.*pi/T*t);
00101 }
00102 for (register int i=0; i<nc; i++) {
00103 DataType alpha = std::atan((xc[i](1)-xsc[1])/(xc[i](0)-xsc[0]));
00104 if (xc[i](0)<xsc[0]) alpha += pi;
00105 DataType xp = xsc[0]+rd*std::cos(alpha)-xs[0];
00106 DataType yp = xsc[1]+rd*std::sin(alpha)-xs[1];
00107 DataType rc = std::sqrt(xp*xp+yp*yp);
00108 aux[i*base::NAux()] = rc*vc[0];
00109 aux[i*base::NAux()+1] = rc*vc[1];
00110 }
00111 }
00112 };
00113
00114
00115 class GFMLevelSetSpecific : public GFMLevelSet<DataType,DIM> {
00116 typedef GFMLevelSet<DataType,DIM> base;
00117
00118 public:
00119 typedef base::grid_fct_type grid_fct_type;
00120 typedef base::grid_data_type grid_data_type;
00121
00122 GFMLevelSetSpecific() {}
00123 virtual ~GFMLevelSetSpecific() {}
00124
00125 virtual void SetGrid(grid_data_type& gdphi, const int& Level, const double& t) {
00126 int Nx = gdphi.extents(0), Ny = gdphi.extents(1);
00127 DCoords lb = base::GH().worldCoords(gdphi.lower(), gdphi.stepsize());
00128 DCoords dx = base::GH().worldStep(gdphi.stepsize());
00129 DataType *phi = (DataType *)gdphi.databuffer();
00130
00131 DataType xsc[2];
00132 xsc[0] = xs[0]; xsc[1] = xs[1];
00133 if (T!=0.) {
00134 xsc[0] += rs*std::cos(2.*pi/T*t);
00135 xsc[1] += rs*std::sin(2.*pi/T*t);
00136 }
00137 for (register int j=0; j<Ny; j++) {
00138 DataType yd = (DataType(j)+static_cast<DataType>(0.5))*dx(1) + lb(1) - xsc[1];
00139 for (register int i=0; i<Nx; i++) {
00140 DataType xd = (DataType(i)+static_cast<DataType>(0.5))*dx(0) + lb(0) - xsc[0];
00141 phi[i+Nx*j] = std::sqrt(xd*xd+yd*yd) - rd;
00142 }
00143 }
00144 }
00145 };
00146
00147
00148 class GhostFluidMethodSpecific : public GhostFluidMethod<VectorType,DIM> {
00149 typedef GhostFluidMethod<VectorType,DIM> base;
00150 public:
00151 GhostFluidMethodSpecific(boundary_type* boundary, levelset_type* levelset) :
00152 base(boundary,levelset) {
00153 xs[0]=0.; xs[1]=0.; rs = 1.; rd=1.; T=0.;
00154 }
00155
00156 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00157 base::register_at(Ctrl,prefix);
00158 RegisterAt(base::LocCtrl,"Center(1)",xs[0]);
00159 RegisterAt(base::LocCtrl,"Center(2)",xs[1]);
00160 RegisterAt(base::LocCtrl,"CircleDuration",T);
00161 RegisterAt(base::LocCtrl,"RadiusCircle",rs);
00162 RegisterAt(base::LocCtrl,"RadiusSphere",rd);
00163 }
00164 };
00165
00166
00167 class SolverSpecific :
00168 public AMRGFMSolver<VectorType,FixupType,FlagType,DIM> {
00169 typedef AMRGFMSolver<VectorType,FixupType,FlagType,DIM> base;
00170 typedef VectorType::InternalDataType DataType;
00171 typedef SchemeFileOutput<SchemeType,DIM> output_type;
00172 public:
00173 SolverSpecific(IntegratorSpecific& integ,
00174 base::initial_condition_type& init,
00175 base::boundary_conditions_type& bc) : base(integ, init, bc) {
00176 SetLevelTransfer(new F77LevelTransfer<VectorType,DIM>(f_prolong, f_restrict));
00177 SetFileOutput(new SchemeGFMFileOutput<SchemeType,FixupType,FlagType,DIM>(*this,integ.Scheme()));
00178 SetFixup(new FixupSpecific());
00179 SetFlagging(new FlaggingSpecific(*this));
00180 AddGFM(new GhostFluidMethodSpecific(new GFMBoundarySpecific(integ.Scheme()),
00181 new GFMLevelSetSpecific()));
00182 }
00183
00184 ~SolverSpecific() {
00185 DeleteGFM(_GFM[0]);
00186 delete _LevelTransfer;
00187 delete _Flagging;
00188 delete _Fixup;
00189 delete _FileOutput;
00190 }
00191 };