00001 #ifdef CH_LANG_CC
00002
00003
00004
00005
00006
00007
00008
00009 #endif
00010
00011 #ifndef _VISCOUSTENSOROP_H_
00012 #define _VISCOUSTENSOROP_H_
00013
00014 #include "AMRMultiGrid.H"
00015 #include "REAL.H"
00016 #include "Box.H"
00017 #include "LevelDataOps.H"
00018 #include "BCFunc.H"
00019 #include "FArrayBox.H"
00020 #include "FluxBox.H"
00021 #include "CFIVS.H"
00022 #include "AMRTGA.H"
00023 #include "TensorCFInterp.H"
00024 #include "CoarseAverage.H"
00025 #include "LevelFluxRegister.H"
00026 #include "NamespaceHeader.H"
00027
00029
00033 class ViscousTensorOp : public LevelTGAHelmOp<LevelData<FArrayBox>, FluxBox>
00034 {
00035 public:
00036
00037 virtual void setAlphaAndBeta(const Real& a_alpha,
00038 const Real& a_beta)
00039 {
00040 m_alpha = a_alpha;
00041 m_beta = a_beta;
00042 defineRelCoef();
00043 }
00044
00045 virtual void diagonalScale(LevelData<FArrayBox>& a_rhs)
00046 {
00047 DisjointBoxLayout grids = a_rhs.disjointBoxLayout();
00048 for(DataIterator dit = grids.dataIterator(); dit.ok(); ++dit)
00049 {
00050 for(int idir = 0; idir < SpaceDim; idir++)
00051 {
00052 int isrc = 0; int idst = idir; int inco = 1;
00053 a_rhs[dit()].mult((*m_acoef)[dit()], isrc, idst, inco);
00054 }
00055 }
00056 }
00057
00059
00061 virtual ~ViscousTensorOp() { ; }
00062
00063
00065
00067 ViscousTensorOp(const DisjointBoxLayout& a_grids,
00068 const DisjointBoxLayout& a_gridsFine,
00069 const DisjointBoxLayout& a_gridsCoar,
00070 const RefCountedPtr<LevelData<FluxBox> >& a_eta,
00071 const RefCountedPtr<LevelData<FluxBox> >& a_lambda,
00072 const RefCountedPtr<LevelData<FArrayBox> >& a_acoef,
00073 Real a_alpha,
00074 Real a_beta,
00075 int a_refToFine,
00076 int a_refToCoar,
00077 const ProblemDomain& a_domain,
00078 const Real& a_dxLevel,
00079 const Real& a_dxCoar,
00080 BCFunc a_bc);
00081
00082 virtual void residual( LevelData<FArrayBox>& a_lhs,
00083 const LevelData<FArrayBox>& a_phi,
00084 const LevelData<FArrayBox>& a_rhs,
00085 bool a_homogeneous = false);
00086
00087 virtual void preCond( LevelData<FArrayBox>& a_correction,
00088 const LevelData<FArrayBox>& a_residual);
00089
00090 virtual void applyOp( LevelData<FArrayBox>& a_lhs,
00091 const LevelData<FArrayBox>& a_phi,
00092 bool a_homogeneous = false);
00093 virtual void create( LevelData<FArrayBox>& a_lhs,
00094 const LevelData<FArrayBox>& a_rhs);
00095 virtual void createCoarsened( LevelData<FArrayBox>& a_lhs,
00096 const LevelData<FArrayBox>& a_rhs,
00097 const int& a_refRat);
00098
00099 void reflux(const LevelData<FArrayBox>& a_phiFine,
00100 const LevelData<FArrayBox>& a_phi,
00101 LevelData<FArrayBox>& residual,
00102 AMRLevelOp<LevelData<FArrayBox> >* a_finerOp);
00103
00104
00105 virtual void getFlux(FluxBox& a_flux,
00106 const LevelData<FArrayBox>& a_data,
00107 const Box& a_grid,
00108 const DataIndex& a_dit,
00109 Real a_scale)
00110 {
00111 const FArrayBox& data = a_data[a_dit];
00112 a_flux.define(a_grid, SpaceDim);
00113 for(int idir = 0; idir < SpaceDim; idir++)
00114 {
00115 const FArrayBox& etaFace = (*m_eta)[a_dit][idir];
00116 const FArrayBox& lambdaFace = (*m_lambda)[a_dit][idir];
00117 const FArrayBox& gradData = m_grad[a_dit];
00118 Box faceBox = a_flux[idir].box();
00119 Box domFaceBox = surroundingNodes(m_domain.domainBox(), idir);
00120 faceBox &= domFaceBox;
00121 getFlux(a_flux[idir], data, gradData, etaFace, lambdaFace, faceBox, idir, 1);
00122 a_flux[idir] *= a_scale;
00123 }
00124 }
00125 void applyOpNoBoundary( LevelData<FArrayBox>& a_lhs,
00126 const LevelData<FArrayBox>& a_phi)
00127 {
00128 this->fillGrad(a_phi);
00129 computeOperatorNoBCs(a_lhs, a_phi);
00130 }
00131
00132 void getFlux(FArrayBox& a_flux,
00133 const FArrayBox& a_data,
00134 const FArrayBox& a_gradData,
00135 const FArrayBox& a_etaFace,
00136 const FArrayBox& a_lambdaFace,
00137 const Box& a_facebox,
00138 int a_dir,
00139 int ref = 1);
00140
00142 void computeOperatorNoBCs(LevelData<FArrayBox>& a_lhs,
00143 const LevelData<FArrayBox>& a_phi);
00144
00145
00146 virtual void assign( LevelData<FArrayBox>& a_lhs,
00147 const LevelData<FArrayBox>& a_rhs) ;
00148 virtual Real dotProduct(const LevelData<FArrayBox>& a_1,
00149 const LevelData<FArrayBox>& a_2) ;
00150 virtual void incr( LevelData<FArrayBox>& a_lhs,
00151 const LevelData<FArrayBox>& a_x,
00152 Real a_scale) ;
00153 virtual void axby( LevelData<FArrayBox>& a_lhs,
00154 const LevelData<FArrayBox>& a_x,
00155 const LevelData<FArrayBox>& a_y,
00156 Real a, Real b) ;
00157
00158 virtual void scale(LevelData<FArrayBox>& a_lhs, const Real& a_scale) ;
00159
00160 virtual Real norm(const LevelData<FArrayBox>& a_x, int a_ord);
00161
00162 virtual void setToZero( LevelData<FArrayBox>& a_x);
00168
00169 virtual void relax(LevelData<FArrayBox>& a_e,
00170 const LevelData<FArrayBox>& a_residual,
00171 int iterations);
00172
00173 virtual void createCoarser(LevelData<FArrayBox>& a_coarse,
00174 const LevelData<FArrayBox>& a_fine,
00175 bool ghosted);
00180 virtual void restrictResidual(LevelData<FArrayBox>& a_resCoarse,
00181 LevelData<FArrayBox>& a_phiFine,
00182 const LevelData<FArrayBox>& a_rhsFine);
00183
00188 virtual void prolongIncrement(LevelData<FArrayBox>& a_phiThisLevel,
00189 const LevelData<FArrayBox>& a_correctCoarse);
00190
00196
00198 virtual int refToCoarser(){ return m_refToCoar; }
00199
00201 virtual void AMRResidual(LevelData<FArrayBox>& a_residual,
00202 const LevelData<FArrayBox>& a_phiFine,
00203 const LevelData<FArrayBox>& a_phi,
00204 const LevelData<FArrayBox>& a_phiCoarse,
00205 const LevelData<FArrayBox>& a_rhs,
00206 bool a_homogeneousPhysBC,
00207 AMRLevelOp<LevelData<FArrayBox> >* a_finerOp);
00208
00211 virtual void AMRResidualNC(LevelData<FArrayBox>& a_residual,
00212 const LevelData<FArrayBox>& a_phiFine,
00213 const LevelData<FArrayBox>& a_phi,
00214 const LevelData<FArrayBox>& a_rhs,
00215 bool a_homogeneousPhysBC,
00216 AMRLevelOp<LevelData<FArrayBox> >* a_finerOp);
00217
00219 virtual void AMRResidualNF(LevelData<FArrayBox>& a_residual,
00220 const LevelData<FArrayBox>& a_phi,
00221 const LevelData<FArrayBox>& a_phiCoarse,
00222 const LevelData<FArrayBox>& a_rhs,
00223 bool a_homogeneousPhysBC);
00224
00231 virtual void AMRRestrict(LevelData<FArrayBox>& a_resCoarse,
00232 const LevelData<FArrayBox>& a_residual,
00233 const LevelData<FArrayBox>& a_correction,
00234 const LevelData<FArrayBox>& a_coarseCorrection);
00235
00237 virtual void AMRProlong(LevelData<FArrayBox>& a_correction,
00238 const LevelData<FArrayBox>& a_coarseCorrection);
00239
00241 virtual void AMRUpdateResidual(LevelData<FArrayBox>& a_residual,
00242 const LevelData<FArrayBox>& a_correction,
00243 const LevelData<FArrayBox>& a_coarseCorrection);
00244
00246
00249 virtual Real AMRNorm(const LevelData<FArrayBox>& a_coarseResid,
00250 const LevelData<FArrayBox>& a_fineResid,
00251 const int& a_refRat,
00252 const int& a_ord);
00253
00254 void homogeneousCFInterp(LevelData<FArrayBox>& a_phif);
00256
00259 void
00260 homogeneousCFInterpPhi(LevelData<FArrayBox>& a_phif,
00261 const DataIndex& a_datInd,
00262 int a_idir,
00263 Side::LoHiSide a_hiorlo);
00264
00266
00270 void homogeneousCFInterpTanGrad(LevelData<FArrayBox>& a_tanGrad,
00271 const LevelData<FArrayBox>& a_phi,
00272 const DataIndex& a_datInd,
00273 int a_idir,
00274 Side::LoHiSide a_hiorlo);
00275
00276 void interpOnIVSHomo(LevelData<FArrayBox>& a_phif,
00277 const DataIndex& a_datInd,
00278 const int a_idir,
00279 const Side::LoHiSide a_hiorlo,
00280 const IntVectSet& a_interpIVS);
00282
00285 void AMROperator( LevelData<FArrayBox>& a_LofPhi,
00286 const LevelData<FArrayBox>& a_phiFine,
00287 const LevelData<FArrayBox>& a_phi,
00288 const LevelData<FArrayBox>& a_phiCoarse,
00289 bool a_homogeneousDomBC,
00290 AMRLevelOp<LevelData<FArrayBox> >* a_finerOp);
00291
00292 void AMROperatorNF( LevelData<FArrayBox>& a_LofPhi,
00293 const LevelData<FArrayBox>& a_phi,
00294 const LevelData<FArrayBox>& a_phiCoarse,
00295 bool a_homogeneousBC);
00296
00297 virtual void AMROperatorNC( LevelData<FArrayBox>& a_LofPhi,
00298 const LevelData<FArrayBox>& a_phiFine,
00299 const LevelData<FArrayBox>& a_phi,
00300 bool a_homogeneousBC,
00301 AMRLevelOp<LevelData<FArrayBox> >* a_finerOp);
00302
00303 void cellGrad(FArrayBox& a_gradPhi,
00304 const FArrayBox& a_phi,
00305 const Box& a_grid);
00306
00307 void cfinterp(const LevelData<FArrayBox>& a_phiFine,
00308 const LevelData<FArrayBox>& a_phiCoarse);
00309
00310 void fillGrad(const LevelData<FArrayBox>& a_phiFine);
00311
00312 static void loHiCenterFace(Box& a_loBox,
00313 int& a_hasLo,
00314 Box& a_hiBox,
00315 int& a_hasHi,
00316 Box& a_centerBox,
00317 const ProblemDomain& a_eblg,
00318 const Box& a_inBox,
00319 const int& a_dir);
00320
00321 static void getFaceDivAndGrad(FArrayBox& a_faceDiv,
00322 FArrayBox& a_faceGrad,
00323 const FArrayBox& a_data,
00324 const FArrayBox& a_gradData,
00325 const ProblemDomain& a_domain,
00326 const Box& a_faceBox,
00327 const int& a_faceDir,
00328 const Real a_dx);
00329
00330 static void getFluxFromDivAndGrad(FArrayBox& a_flux,
00331 const FArrayBox& a_faceDiv,
00332 const FArrayBox& a_faceGrad,
00333 const FArrayBox& a_etaFace,
00334 const FArrayBox& a_lambdaFace,
00335 const Box& a_faceBox,
00336 int a_dir);
00337
00339
00344 void divergenceCC(LevelData<FArrayBox>& a_div,
00345 const LevelData<FArrayBox>& a_phi,
00346 const LevelData<FArrayBox>* a_phiC);
00347 protected:
00348 void defineRelCoef();
00349 RefCountedPtr<LevelData<FluxBox> > m_eta;
00350 RefCountedPtr<LevelData<FluxBox> > m_lambda;
00351 RefCountedPtr<LevelData<FArrayBox> > m_acoef;
00352 Real m_alpha;
00353 Real m_beta;
00354 int m_refToCoar;
00355 int m_refToFine;
00356 BCFunc m_bc;
00357
00358
00359 int m_ncomp;
00360 Real m_dx;
00361 Real m_dxCrse;
00362 ProblemDomain m_domain;
00363
00364 LevelData<FArrayBox> m_relaxCoef;
00365 LevelData<FArrayBox> m_grad;
00366 LevelDataOps<FArrayBox> m_levelOps;
00367 Copier m_exchangeCopier;
00368 TensorCFInterp m_interpWithCoarser;
00369
00370 LayoutData<CFIVS> m_loCFIVS[SpaceDim];
00371 LayoutData<CFIVS> m_hiCFIVS[SpaceDim];
00372
00373 LayoutData<TensorFineStencilSet> m_hiTanStencilSets[SpaceDim];
00374 LayoutData<TensorFineStencilSet> m_loTanStencilSets[SpaceDim];
00375 Vector<IntVect> m_colors;
00376
00377 private:
00379 ViscousTensorOp()
00380 {
00381 MayDay::Error("invalid operator");
00382 }
00383
00384 ViscousTensorOp(const ViscousTensorOp& a_opin)
00385 {
00386 MayDay::Error("invalid operator");
00387 }
00388 void operator=(const ViscousTensorOp& a_opin)
00389 {
00390 MayDay::Error("invalid operator");
00391 }
00392 };
00393
00395
00398 class ViscousTensorOpFactory: public AMRLevelOpFactory<LevelData<FArrayBox> >
00399 {
00400 public:
00401 virtual ~ViscousTensorOpFactory(){;}
00402
00403 ViscousTensorOpFactory(const Vector<DisjointBoxLayout>& a_grids,
00404 const Vector<RefCountedPtr<LevelData<FluxBox> > >& a_eta,
00405 const Vector<RefCountedPtr<LevelData<FluxBox> > >& a_lambda,
00406 const Vector<RefCountedPtr<LevelData<FArrayBox> > >& a_acoef,
00407 Real a_alpha,
00408 Real a_beta,
00409 const Vector<int>& a_refRatios,
00410 const ProblemDomain& a_domainCoar,
00411 const Real& a_dxCoar,
00412 BCFunc a_bc);
00413
00415 virtual ViscousTensorOp*
00416 MGnewOp(const ProblemDomain& a_FineindexSpace,
00417 int depth,
00418 bool homoOnly = true);
00419
00421 virtual ViscousTensorOp* AMRnewOp(const ProblemDomain& a_indexSpace);
00422
00424 virtual int refToFiner(const ProblemDomain&) const;
00425
00426 private:
00427 Vector<RefCountedPtr<LevelData<FluxBox> > > m_eta;
00428 Vector<RefCountedPtr<LevelData<FluxBox> > > m_lambda;
00429 Vector<RefCountedPtr<LevelData<FArrayBox> > > m_acoef;
00430 LevelData<FArrayBox> m_phic;
00431 Vector<ProblemDomain> m_domains;
00432 Vector<DisjointBoxLayout> m_boxes;
00433 Vector<Real> m_dx;
00434
00435 Vector<int> m_refRatios;
00436 BCFunc m_bc;
00437 Real m_alpha;
00438 Real m_beta;
00439
00441 ViscousTensorOpFactory()
00442 {
00443 MayDay::Error("invalid operator");
00444 }
00445
00446 ViscousTensorOpFactory(const ViscousTensorOpFactory& a_opin)
00447 {
00448 MayDay::Error("invalid operator");
00449 }
00450
00451 void operator=(const ViscousTensorOpFactory& a_opin)
00452 {
00453 MayDay::Error("invalid operator");
00454 }
00455 };
00456
00457 extern void
00458 coarsenStuff(LevelData<FluxBox> & a_etaCoar,
00459 LevelData<FluxBox> & a_lambdaCoar,
00460 const LevelData<FluxBox> & a_etaFine,
00461 const LevelData<FluxBox> & a_lambdaFine,
00462 const int & a_refToDepth);
00463
00464 #include "NamespaceFooter.H"
00465 #endif
