vtf-logo

clawpack/applications/euler/2d/VortexRotGFM/src/Problem.h

00001 // -*- C++ -*-
00002 
00003 // Copyright (C) 2003-2007 California Institute of Technology
00004 // Ralf Deiterding, ralf@amroc.net
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     // Append to output file only on processor VizServer
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;