releasedExamples/MultiDimPhase/advectSrc/OldPatchGodunov.cpp
// Compute the flattening coefficients from the primitive variables
Compute the flattening coefficients from the primitive variables
Compute the flattening coefficients from the primitive variables
void OldPatchGodunov::computeFlattening(FArrayBox&       a_flattening,
                                        const FArrayBox& a_W,
                                        const Box&       a_box)
{
  CH_assert(isDefined());
  CH_assert(a_W.box().contains(a_box));
  int idir;
  FArrayBox zetaDir(a_box,SpaceDim);
  FArrayBox dVel(a_box,SpaceDim);
  Interval velInterval= velocityInterval();
  int v0index = velInterval.begin();
  for (idir = 0; idir < SpaceDim; idir++)
  {
    Box box1 = a_box;
    box1.grow(idir,1);
    Box box2 = a_box;
    box2.grow(idir,2);
    Box box3 = a_box;
    box3.grow(idir,3);
    Box loBox,hiBox,centerBox,entireBox;
    int hasLo,hasHi;
    oldLoHiCenter(loBox,hasLo,hiBox,hasHi,centerBox,entireBox,
                  box3,m_domain,idir);
    FArrayBox delta1p(entireBox,1);
    int pressIndex = pressureIndex();
    FORT_OLDGETGRADF(CHF_FRA1(delta1p,0),
                     CHF_CONST_FRA1(a_W,pressIndex),
                     CHF_CONST_INT(idir),
                     CHF_BOX(loBox),
                     CHF_CONST_INT(hasLo),
                     CHF_BOX(hiBox),
                     CHF_CONST_INT(hasHi),
                     CHF_BOX(centerBox));
    oldLoHiCenter(loBox,hasLo,hiBox,hasHi,centerBox,entireBox,
                  box2,m_domain,idir);
    FArrayBox delta2p(entireBox,1);
    FORT_OLDGETDPTWOF(CHF_FRA1(delta2p,0),
                      CHF_CONST_FRA1(delta1p,0),
                      CHF_CONST_INT(idir),
                      CHF_BOX(loBox),
                      CHF_CONST_INT(hasLo),
                      CHF_BOX(hiBox),
                      CHF_CONST_INT(hasHi),
                      CHF_BOX(centerBox));
    FArrayBox bulkMin(entireBox,1);
    int bulkIndex = bulkModulusIndex();
    FORT_OLDMIN3PTSF(CHF_FRA1(bulkMin,0),
                     CHF_CONST_FRA1(a_W,bulkIndex),
                     CHF_CONST_INT(idir),
                     CHF_BOX(loBox),
                     CHF_CONST_INT(hasLo),
                     CHF_BOX(hiBox),
                     CHF_CONST_INT(hasHi),
                     CHF_BOX(centerBox));
    FArrayBox zetaTwiddleDir(entireBox,1);
    FORT_OLDGETFLATF(CHF_FRA1(zetaTwiddleDir,0),
                     CHF_CONST_FRA1(delta1p,0),
                     CHF_CONST_FRA1(delta2p,0),
                     CHF_CONST_FRA1(bulkMin,0),
                     CHF_BOX(entireBox));
    oldLoHiCenter(loBox,hasLo,hiBox,hasHi,centerBox,entireBox,
                  box1,m_domain,idir);
    FORT_OLDMIN3PTSF(CHF_FRA1(zetaDir,idir),
                     CHF_CONST_FRA1(zetaTwiddleDir,0),
                     CHF_CONST_INT(idir),
                     CHF_BOX(loBox),
                     CHF_CONST_INT(hasLo),
                     CHF_BOX(hiBox),
                     CHF_CONST_INT(hasHi),
                     CHF_BOX(centerBox));
    FORT_OLDGETGRADF(CHF_FRA1(dVel,idir),
                     CHF_CONST_FRA1(a_W,v0index+idir),
                     CHF_CONST_INT(idir),
                     CHF_BOX(loBox),
                     CHF_CONST_INT(hasLo),
                     CHF_BOX(hiBox),
                     CHF_CONST_INT(hasHi),
                     CHF_BOX(centerBox));
  }
  FORT_OLDMINFLATF(CHF_FRA1(a_flattening,0),
                   CHF_CONST_FRA(zetaDir),
                   CHF_CONST_FRA(dVel),
                   CHF_BOX(a_box));
}
void OldPatchGodunov :: computeFlattening ( FArrayBox & a_flattening , const FArrayBox & a_W , const Box & a_box ) { CH_assert ( isDefined ( ) ) ; CH_assert ( a_W . box ( ) . contains ( a_box ) ) ; int idir ; FArrayBox zetaDir ( a_box , SpaceDim ) ; FArrayBox dVel ( a_box , SpaceDim ) ; Interval velInterval = velocityInterval ( ) ; int v0index = velInterval . begin ( ) ; for ( idir = 0 ; idir < SpaceDim ; idir ++ ) { Box box1 = a_box ; box1 . grow ( idir , 1 ) ; Box box2 = a_box ; box2 . grow ( idir , 2 ) ; Box box3 = a_box ; box3 . grow ( idir , 3 ) ; Box loBox , hiBox , centerBox , entireBox ; int hasLo , hasHi ; oldLoHiCenter ( loBox , hasLo , hiBox , hasHi , centerBox , entireBox , box3 , m_domain , idir ) ; FArrayBox delta1p ( entireBox , 1 ) ; int pressIndex = pressureIndex ( ) ; FORT_OLDGETGRADF ( CHF_FRA1 ( delta1p , 0 ) , CHF_CONST_FRA1 ( a_W , pressIndex ) , CHF_CONST_INT ( idir ) , CHF_BOX ( loBox ) , CHF_CONST_INT ( hasLo ) , CHF_BOX ( hiBox ) , CHF_CONST_INT ( hasHi ) , CHF_BOX ( centerBox ) ) ; oldLoHiCenter ( loBox , hasLo , hiBox , hasHi , centerBox , entireBox , box2 , m_domain , idir ) ; FArrayBox delta2p ( entireBox , 1 ) ; FORT_OLDGETDPTWOF ( CHF_FRA1 ( delta2p , 0 ) , CHF_CONST_FRA1 ( delta1p , 0 ) , CHF_CONST_INT ( idir ) , CHF_BOX ( loBox ) , CHF_CONST_INT ( hasLo ) , CHF_BOX ( hiBox ) , CHF_CONST_INT ( hasHi ) , CHF_BOX ( centerBox ) ) ; FArrayBox bulkMin ( entireBox , 1 ) ; int bulkIndex = bulkModulusIndex ( ) ; FORT_OLDMIN3PTSF ( CHF_FRA1 ( bulkMin , 0 ) , CHF_CONST_FRA1 ( a_W , bulkIndex ) , CHF_CONST_INT ( idir ) , CHF_BOX ( loBox ) , CHF_CONST_INT ( hasLo ) , CHF_BOX ( hiBox ) , CHF_CONST_INT ( hasHi ) , CHF_BOX ( centerBox ) ) ; FArrayBox zetaTwiddleDir ( entireBox , 1 ) ; FORT_OLDGETFLATF ( CHF_FRA1 ( zetaTwiddleDir , 0 ) , CHF_CONST_FRA1 ( delta1p , 0 ) , CHF_CONST_FRA1 ( delta2p , 0 ) , CHF_CONST_FRA1 ( bulkMin , 0 ) , CHF_BOX ( entireBox ) ) ; oldLoHiCenter ( loBox , hasLo , hiBox , hasHi , centerBox , entireBox , box1 , m_domain , idir ) ; FORT_OLDMIN3PTSF ( CHF_FRA1 ( zetaDir , idir ) , CHF_CONST_FRA1 ( zetaTwiddleDir , 0 ) , CHF_CONST_INT ( idir ) , CHF_BOX ( loBox ) , CHF_CONST_INT ( hasLo ) , CHF_BOX ( hiBox ) , CHF_CONST_INT ( hasHi ) , CHF_BOX ( centerBox ) ) ; FORT_OLDGETGRADF ( CHF_FRA1 ( dVel , idir ) , CHF_CONST_FRA1 ( a_W , v0index + idir ) , CHF_CONST_INT ( idir ) , CHF_BOX ( loBox ) , CHF_CONST_INT ( hasLo ) , CHF_BOX ( hiBox ) , CHF_CONST_INT ( hasHi ) , CHF_BOX ( centerBox ) ) ; } FORT_OLDMINFLATF ( CHF_FRA1 ( a_flattening , 0 ) , CHF_CONST_FRA ( zetaDir ) , CHF_CONST_FRA ( dVel ) , CHF_BOX ( a_box ) ) ; }
Compute the flattening coefficients from the primitive variables
Compute the flattening coefficients from the primitive variables
// The current direction // The directional flattening coefficients // The divergence of the velocity // The interval of the primitive variables corresponding to the velocity // Get the directional flattening coefficients in each direction // A box one larger (in direction "idir") than the final result box // A box two larger (in direction "idir") than the final result box // A box three larger (in direction "idir") than the final result box // Compute where centered differences can be used and where one sided // differences need to be used.  The data used for the differences is // defined on "box3" // Compute the first differences in "pressure" // Compute where centered differences can be used and where one sided // differences need to be used.  The data used for the differences is // defined on "box2" // Compute the second differences in "pressure" // Compute a 3-way minimum of the "bulk modulus" // Use the first and second differences normalized by the 3-way minimum // computed above to generate directional flattening coefficients // Compute where centered differences can be used and where one sided // differences need to be used.  The data used for the differences is // defined on "box1" // Take a 3-way minimum of the directional flattening coefficients // Compute each component of the divergence of the velocity // At each point, set the flattening coefficient to the minimum of all // the directional flattening coefficients if the divergence of the velocity // is negative, otherwise set it to 1 (no flattening).