00001
00002
00003
00004
00005
00006 #ifndef AMROC_PROBLEM_H
00007 #define AMROC_PROBLEM_H
00008
00009 #include "euler2.h"
00010 #include "ClpProblem.h"
00011
00012 #define f_ipvel FORTRAN_NAME(ipvel, IPVEL)
00013
00014 extern "C" {
00015 void f_ipvel();
00016 }
00017
00018 #define OWN_GFMAMRSOLVER
00019 #include "ClpStdGFMProblem.h"
00020 #include "F77Interfaces/F77GFMExactSolution.h"
00021
00022 class SolverSpecific :
00023 public AMRGFMSolver<VectorType,FixupType,FlagType,DIM> {
00024 typedef VectorType::InternalDataType DataType;
00025 typedef AMRGFMSolver<VectorType,FixupType,FlagType,DIM> base;
00026 typedef GhostFluidMethod<VectorType,DIM>::boundary_conditions_type gfm_bc_type;
00027 public:
00028 SolverSpecific(IntegratorSpecific& integ,
00029 base::initial_condition_type& init,
00030 base::boundary_conditions_type& bc) : base(integ, init, bc) {
00031 SetLevelTransfer(new F77LevelTransfer<VectorType,DIM>(f_prolong, f_restrict));
00032 #ifdef f_flgout
00033 SetFileOutput(new F77GFMFileOutput<VectorType,FixupType,FlagType,DIM>(*this,f_flgout));
00034 #else
00035 SetFileOutput(new GFMFileOutput<VectorType,FixupType,FlagType,DIM>(*this));
00036 #endif
00037 SetFixup(new FixupSpecific(integ));
00038 SetFlagging(new FlaggingSpecific(*this));
00039 AddGFM(new GhostFluidMethod<VectorType,DIM>(
00040 new F77GFMBoundary<VectorType,DIM>(f_ibndrfl,f_itrans,f_ipvel),
00041 new F77GFMLevelSet<DataType,DIM>(f_lset),
00042 (gfm_bc_type*)(new F77BoundaryConditions<Vector<DataType,1>,DIM>(f_boundary))));
00043 SetExactSolution(new F77GFMExactSolution<VectorType,FixupType,FlagType,DIM>(*this,f_exact));
00044 }
00045
00046 ~SolverSpecific() {
00047 delete _GFM[0]->GetBoundaryConditionsP();
00048 DeleteGFM(_GFM[0]);
00049 delete _Flagging;
00050 delete _Fixup;
00051 }
00052
00053 virtual void SetupData() {
00054 base::SetupData();
00055 _GFM[0]->BoundaryConditions_().SetupData(base::PGH(), base::NGhosts());
00056 }
00057
00058 virtual void Initialize_(const double& dt_start) {
00059 base::Initialize_(dt_start);
00060 CalculateCurrentMass(mass_old);
00061 }
00062
00063 virtual double Tick(int VariableTimeStepping, const double dtv[],
00064 const double cflv[], int& Rejections) {
00065
00066 double cfl = base::Tick(VariableTimeStepping,dtv,cflv,Rejections);
00067
00068 VectorType mass_new;
00069 CalculateCurrentMass(mass_new);
00070
00071
00072 int me = MY_PROC;
00073 if (me == VizServer) {
00074 std::ofstream outfile;
00075 std::ostream* out;
00076 std::string fname = "mass.dat";
00077 outfile.open(fname.c_str(), std::ios::out | std::ios::app);
00078 out = new std::ostream(outfile.rdbuf());
00079 *out << base::t[0];
00080 for (int n=0; n<base::NEquations(); n++)
00081 *out << " " << mass_old[n]-mass_new[n];
00082 *out << std::endl;
00083 outfile.close();
00084 delete out;
00085 }
00086
00087 return cfl;
00088 }
00089
00090 void CalculateCurrentMass(VectorType& mass) {
00091 int Level=0;
00092 int Time = CurrentTime(base::GH(),Level);
00093 int TStep = TimeStep(base::U(),Level);
00094 forall (base::U(),Time,Level,c)
00095 base::U()(Time+TStep,Level,c).equals(base::U()(Time,Level,c));
00096 end_forall
00097
00098 if (base::NGFM()>0) {
00099 forall (base::U(),Time+TStep,Level,c)
00100 for (int g=0; g<base::NGFM(); g++) {
00101 if (!base::GFM(g).IsUsed()) continue;
00102 BeginFastIndex2(phi, base::GFM(g).Phi()(Time,Level,c).bbox(),
00103 base::GFM(g).Phi()(Time,Level,c).data(), DataType);
00104 BeginFastIndex2(u, base::U()(Time+TStep,Level,c).bbox(),
00105 base::U()(Time+TStep,Level,c).data(), VectorType);
00106 for_2 (n, m, base::GFM(g).Phi()(Time,Level,c).bbox(),
00107 base::GFM(g).Phi()(Time,Level,c).bbox().stepsize())
00108 if (FastIndex2(phi,n,m)<-base::GFM(g).Boundary().Cutoff())
00109 FastIndex2(u,n,m) = static_cast<VectorType>(0.0);
00110 end_for
00111 EndFastIndex2(u);
00112 EndFastIndex2(phi);
00113 }
00114 end_forall
00115 }
00116
00117 DCoords dx = base::GH().worldStep(base::U().GF_StepSize(Level));
00118 mass = Sum(base::U(),Time+TStep,Level);
00119 for (int d=0; d<base::Dim(); d++) mass *= dx(d);
00120 }
00121
00122 protected:
00123 VectorType mass_old;