00001
00002
00003
00004
00005
00006 #ifndef AMROCFLUIDPROBLEMH
00007 #define AMROCFLUIDPROBLEMH
00008
00009 #include "eulerznd3.h"
00010
00011 #undef NAUX
00012 #define NAUX 2
00013
00014 #include "ClpProblem.h"
00015
00016 #define MaxIntPoints (20)
00017 #define MaxDPoints (10)
00018
00019 #define OWNFLAGGING
00020 #define OWNGFMAMRSOLVER
00021 #include "ClpStdGFMProblem.h"
00022 #include "AMRGFMInterpolation.h"
00023 #include "StatCPTLevelSet.h"
00024
00025 class FlaggingSpecific :
00026 public AMRFlagging<VectorType,FixupType,FlagType,DIM> {
00027 typedef AMRFlagging<VectorType,FixupType,FlagType,DIM> base;
00028 public:
00029 FlaggingSpecific(base::solvertype& solver) : base(solver) {
00030 base::AddCriterion(new ByValue<VectorType,FlagType,DIM>());
00031 base::AddCriterion(new ScaledGradient<VectorType,FlagType,DIM>());
00032 base::AddCriterion(new LimiterType<VectorType,FlagType,DIM>());
00033 base::AddCriterion(new AbsoluteError<VectorType,FixupType,FlagType,DIM>(solver));
00034 base::AddCriterion(new RelativeError<VectorType,FixupType,FlagType,DIM>(solver));
00035 #ifdef fflgout
00036 base::AddCriterion(new F77ByValue<VectorType,FlagType,DIM>(fflgout));
00037 base::AddCriterion(new F77ScaledGradient<VectorType,FlagType,DIM>(fflgout));
00038 base::AddCriterion(new F77LimiterType<VectorType,FlagType,DIM>(fflgout));
00039 base::AddCriterion(new F77AbsoluteError<VectorType,FixupType,FlagType,DIM>(solver,fflgout));
00040 base::AddCriterion(new F77RelativeError<VectorType,FixupType,FlagType,DIM>(solver,fflgout));
00041 #endif
00042 dcmin = new int[MaxDPoints*base::Dim()];
00043 dcmax = new int[MaxDPoints*base::Dim()];
00044 dclev = new int[MaxDPoints];
00045 for (int d=0; d<MaxDPoints*base::Dim(); d++) {
00046 dcmin[d] = 0; dcmax[d] = -1;
00047 }
00048 for (int i=0; i<MaxDPoints; i++)
00049 dclev[i] = 1000000;
00050 }
00051
00052 ~FlaggingSpecific() {
00053 DeleteAllCriterions();
00054 delete [] dcmin;
00055 delete [] dcmax;
00056 }
00057
00058 virtual void registerat(ControlDevice& Ctrl, const std::string& prefix) {
00059 base::registerat(Ctrl, prefix);
00060 char VariableName[32];
00061 for (int i=0; i<MaxDPoints; i++) {
00062 for (int d=0; d<base::Dim(); d++) {
00063 sprintf(VariableName,"DCellsMin(%d,%d)",i+1,d+1);
00064 RegisterAt(base::LocCtrl,VariableName,dcmin[i*base::Dim()+d]);
00065 sprintf(VariableName,"DCellsMax(%d,%d)",i+1,d+1);
00066 RegisterAt(base::LocCtrl,VariableName,dcmax[i*base::Dim()+d]);
00067 }
00068 sprintf(VariableName,"DMinLevel(%d)",i+1);
00069 RegisterAt(base::LocCtrl,VariableName,dclev[i]);
00070 }
00071 }
00072 virtual void registerat(ControlDevice& Ctrl) {
00073 registerat(Ctrl, "");
00074 }
00075
00076 virtual void SetFlags(const int Time, const int Level, double t, double dt) {
00077 base::SetFlags(Time, Level, t, dt);
00078 for (int i=0; i<MaxDPoints; i++)
00079 if (Level>=dclev[i]) {
00080 Coords lb(base::Dim(), dcmin+i*base::Dim()), ub(base::Dim(), dcmax+i*base::Dim()),
00081 s(base::Dim(),1);
00082 BBox bbox(lb*RefinedBy(base::GH(),MaxLevel(base::GH())),
00083 ub*RefinedBy(base::GH(),MaxLevel(base::GH())),
00084 s*StepSize(base::GH(),Level));
00085 forall (base::Flags(),Time,Level,c)
00086 base::Flags()(Time,Level,c).equals(GoodPoint, bbox);
00087 endforall
00088 }
00089 }
00090
00091 private:
00092 int *dcmin, *dcmax, *dclev;
00093 };
00094
00095
00096 template <class DataType, int dim>
00097 class F77CPTLevelSet : public StatCPTLevelSet<DataType,dim> {
00098 typedef StatCPTLevelSet<DataType,dim> base;
00099 public:
00100 typedef typename base::gridfcttype gridfcttype;
00101 typedef typename base::griddatatype griddatatype;
00102 typedef genericfortranfunc genericfunctype;
00103
00104 typedef void (*lset3functype) ( const INTEGER& maxmx, const INTEGER& maxmy, const INTEGER& maxmz,
00105 const INTEGER& mbc,
00106 const INTEGER& mx, const INTEGER& my, const INTEGER& mz,
00107 const DOUBLE x[], const DOUBLE y[], const DOUBLE z[],
00108 const DOUBLE& dx, const DOUBLE& dy, const DOUBLE& dz,
00109 DataType q[], const DOUBLE& t);
00110
00111 F77CPTLevelSet() : base(), flset(0) {}
00112 F77CPTLevelSet(genericfunctype lset) : base(), flset(lset) {}
00113
00114 virtual void SetPhi(gridfcttype& phi, const int Time, const int Level, double t) {
00115 base::SetPhi(phi,Time,Level,t);
00116 forall (phi,Time,Level,c)
00117 SetGrid(phi(Time,Level,c),Level,t);
00118 endforall
00119 }
00120
00121 virtual void SetGrid(griddatatype& gdphi, const int& level, const double& t) {
00122 assert(flset != (genericfunctype) 0);
00123 Coords ex = gdphi.extents();
00124 DCoords lbcorner = base::GH().worldCoords(gdphi.lower(), gdphi.stepsize());
00125 DCoords dx = base::GH().worldStep(gdphi.stepsize());
00126 int maxm[3], mx[3], d;
00127 DataType* x[3];
00128 for (d=0; d<dim; d++) {
00129 maxm[d] = ex(d);
00130 mx[d] = ex(d)-2*base::NGhosts();
00131 x[d] = new DataType[maxm[d]];
00132 for (int i=0; i<maxm[d]; i++) x[d][i] = (i+0.5)*dx(d)+lbcorner(d);
00133 }
00134
00135 ((lset3functype) flset)(AA(3,mx), base::NGhosts(), AA(3,mx),
00136 AA(3,x), AA(3,dx), gdphi.data(), t);
00137
00138 for (d=0; d<dim; d++)
00139 delete [] x[d];
00140 }
00141
00142 inline void SetFunc(genericfunctype lset) { flset = lset; }
00143 genericfunctype GetFunc() const { return flset; }
00144
00145 protected:
00146 genericfunctype flset;
00147 };
00148
00149
00150
00151 class SolverSpecific :
00152 public AMRGFMSolver<VectorType,FixupType,FlagType,DIM> {
00153 typedef VectorType::InternalDataType DataType;
00154 typedef AMRGFMSolver<VectorType,FixupType,FlagType,DIM> base;
00155 typedef AMRGFMInterpolation<VectorType,FixupType,FlagType,DIM> interpolationtype;
00156 public:
00157 typedef F77GFMFileOutput<VectorType,FixupType,FlagType,DIM> outputtype;
00158 typedef interpolationtype::pointtype pointtype;
00159
00160 SolverSpecific(IntegratorSpecific& integ,
00161 base::initialconditiontype& init,
00162 base::boundaryconditionstype& bc) : base(integ, init, bc) {
00163 SetLevelTransfer(new F77LevelTransfer<VectorType,DIM>(fprolong, frestrict));
00164 SetFileOutput(new outputtype(*this,fflgout));
00165 SetFixup(new FixupSpecific(IntegratorSpecific));
00166 SetFlagging(new FlaggingSpecific(*this));
00167 AddGFM(new GhostFluidMethod<VectorType,DIM>(
00168 new F77GFMBoundary<VectorType,DIM>(fibndrfl,fitrans),
00169 new F77CPTLevelSet<DataType,DIM>(flset)));
00170 Interpolation = new interpolationtype(*this);
00171 NIntPoints = 0;
00172 IntName = "ptrack.txt";
00173 }
00174
00175 ~SolverSpecific() {
00176 DeleteGFM(GFM[0]);
00177 delete Flagging;
00178 delete Fixup;
00179 delete Interpolation;
00180 }
00181
00182 virtual void registerat(ControlDevice& Ctrl, const std::string& prefix) {
00183 base::registerat(Ctrl,prefix);
00184 IntCtrl = base::LocCtrl.getSubDevice("TrackPressure");
00185 RegisterAt(IntCtrl,"NPoints",NIntPoints);
00186 char VariableName[32];
00187 for (int nc=0; nc<MaxIntPoints; nc++) {
00188 for (int d=0; d<base::Dim(); d++) {
00189 std::sprintf(VariableName,"Point(%d,%d)",nc+1,d+1);
00190 RegisterAt(IntCtrl,VariableName,IntPoints[nc](d));
00191 IntPoints[nc](d) = 0.0;
00192 }
00193 }
00194 RegisterAt(IntCtrl,"FileName",IntName);
00195 }
00196 virtual void registerat(ControlDevice& Ctrl) {
00197 base::registerat(Ctrl);
00198 }
00199
00200 virtual void SetupData() {
00201 base::SetupData();
00202 Interpolation->SetupData(base::PGH(), base::NGhosts());
00203 }
00204
00205 virtual void Advance(double& t, double& dt) {
00206 base::Advance(t,dt);
00207 if (NIntPoints<=0 || NIntPoints>MaxIntPoints) return;
00208
00209 int me = MYPROC;
00210
00211 for (register int l=0; l<=FineLevel(base::GH()); l++) {
00212 int Time = CurrentTime(base::GH(),l);
00213 int presscnt = base::Dim()+4;
00214 ((outputtype*) FileOutput)->Transform(base::U(), base::Work(), Time, l,
00215 presscnt, base::t[l]);
00216 }
00217
00218 DataType* p = new DataType[NIntPoints];
00219 Interpolation->PointsValues(base::Work(),base::t[0],NIntPoints,IntPoints,p);
00220 Interpolation->ArrayCombine(VizServer,NIntPoints,p);
00221
00222
00223 if (me == VizServer) {
00224 std::ofstream outfile;
00225 std::ostream* out;
00226 outfile.open(IntName.cstr(), std::ios::out | std::ios::app);
00227 out = new std::ostream(outfile.rdbuf());
00228 *out << base::t[0];
00229 for (int ns=0; ns<NIntPoints; ns++)
00230 *out << " " << p[ns];
00231 *out << std::endl;
00232 outfile.close();
00233 delete out;
00234 }
00235
00236 delete [] p;
00237 }
00238
00239 protected:
00240 IntegratorSpecific IntegratorSpecific;
00241 InitialConditionSpecific InitialConditionSpecific;
00242 BoundaryConditionsSpecific BoundaryConditionsSpecific;
00243 interpolationtype* Interpolation;
00244 ControlDevice IntCtrl;
00245 int NIntPoints;
00246 pointtype IntPoints[MaxIntPoints];
00247 std::string IntName;