00001
00002
00003
00004
00005
00006 #ifndef AMROC_WENO_FIXUP_H
00007 #define AMROC_WENO_FIXUP_H
00008
00016 #include "AMRFixup.h"
00017
00018 template <class VectorType, class FixupType, int dim>
00019 class WENOFixup : public AMRFixup<VectorType,FixupType,dim> {
00020 typedef typename VectorType::InternalDataType DataType;
00021 typedef AMRFixup<VectorType,FixupType,dim> base;
00022
00023 public:
00024 typedef typename base::vec_grid_data_type vec_grid_data_type;
00025
00026 WENOFixup() : base() {}
00027
00028 virtual void SaveFluxes(const int Time, const int Level, const int c,
00029 vec_grid_data_type* flux[],
00030 const double dt, const int& mpass) {
00031
00032 if (!base::U().has_children(Time,Level,c)) return;
00033 DCoords dx = base::GH().worldStep(base::U().GF_StepSize(Level));
00034 const int refinefactor = RefineFactor(base::GH(),Level);
00035
00036 DataType Fac = -dt;
00037 int d;
00038 for (d=0; d<dim; d++)
00039 Fac *= dx(d);
00040
00041 #ifdef DEBUG_PRINT_WHOLE_FIXUP
00042 ( comm_service::log() << "\n*** Saving Coarse Fluxes - "
00043 << base::U().interiorbbox(Time,Level,c) << " ***\n").flush();
00044 #endif
00045 for (register int j=0; j<base::U().children(Time,Level,c); j++) {
00046 BBox work(base::U().childbox(Time,Level,c,j)*base::U().interiorbbox(Time,Level,c));
00047 if (!work.empty()) {
00048 const int& idx = base::U().childidx(Time,Level,c,j);
00049 for (d=0; d<2*dim; d++)
00050 if (base::ValidSide(Time,Level+1,idx,d)) {
00051 BBox bb(base::FixupBBox(coarsen(base::U().interiorbbox(Time,Level+1,idx),refinefactor),d)*work);
00052 if (!bb.empty()) {
00053 bb = base::FluxBBox(bb,d);
00054 #ifdef DEBUG_PRINT_WHOLE_FIXUP
00055 ( comm_service::log() << "To " << base::U().interiorbbox(Time,Level+1,idx)
00056 << " Side: " << d << " using " << bb << " into "
00057 << base::F(d)(Time,Level+1,idx).bbox() << "\n").flush();
00058 #endif
00059 (*flux[d]).multiply(Fac/dx(d/2), bb);
00060
00061 if (mpass<=1)
00062 base::copyf_to(base::F(d)(Time,Level+1,idx),(*flux[d]),bb,d);
00063 else
00064 base::addf_to(base::F(d)(Time,Level+1,idx),(*flux[d]),bb,d);
00065
00066 (*flux[d]).divide(Fac/dx(d/2), bb);
00067 #ifdef DEBUG_PRINT_WHOLE_FIXUP
00068 AlignBBox(bb,d);
00069 base::debug_print_ld(base::F(d)(Time,Level+1,idx),bb);
00070 ( comm_service::log() << "\n" ).flush();
00071 #endif
00072 }
00073 }
00074 }
00075 }
00076 }
00077
00078 virtual void AddIntercellFluxes(const int Time, const int Level, const int c,
00079 vec_grid_data_type* flux[], const double dt,
00080 const int& mpass) {
00081
00082 DCoords dx = base::GH().worldStep(base::U().GF_StepSize(Level));
00083 DataType Fac = dt;
00084 int d;
00085 for (d=0; d<dim; d++)
00086 Fac *= dx(d);
00087
00088 #ifdef DEBUG_PRINT_WHOLE_FIXUP
00089 ( comm_service::log() << "\n*** Adding Higher Order Fluxes - "
00090 << base::U().interiorbbox(Time,Level,c) << " ***\n").flush();
00091 #endif
00092 for (d=0; d<2*dim; d++)
00093 if (base::ValidSide(Time,Level,c,d)) {
00094 BBox bb(base::FluxBBox(base::U().interiorbbox(Time,Level,c), d));
00095 (*flux[d]).multiply(Fac/dx(d/2), bb);
00096 base::addf_to(base::F(d)(Time,Level,c),(*flux[d]),bb,d);
00097 (*flux[d]).divide(Fac/dx(d/2), bb);
00098 #ifdef DEBUG_PRINT_WHOLE_FIXUP
00099 base::debug_print_ld(base::F(d)(Time,Level,c));
00100 ( comm_service::log() << "\n").flush();
00101 #endif
00102 }
00103 }
00104
00105 };
00106
00107 #endif