Skip to content

Commit

Permalink
fix interp functions for single precison (AMReX-Codes#887)
Browse files Browse the repository at this point in the history
  • Loading branch information
@WeiqunZhang
WeiqunZhang authored May 15, 2020
1 parent e5d0c10 commit a017d59
Show file tree
Hide file tree
Showing 3 changed files with 138 additions and 138 deletions.
50 changes: 25 additions & 25 deletions Src/AmrCore/AMReX_Interp_1D_C.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,8 +33,8 @@ ccinterp_compute_voff (Box const& cbx, IntVect const& ratio, Geometry const& cge
const int ic = amrex::coarsen(i, ratio[0]);
const int ii = i - flo.x;
const int iic = ic - clo.x;
const Real fcen = 0.5*(fvc[ii ]+fvc[ii +1]);
const Real ccen = 0.5*(cvc[iic]+cvc[iic+1]);
const Real fcen = 0.5_rt*(fvc[ii ]+fvc[ii +1]);
const Real ccen = 0.5_rt*(cvc[iic]+cvc[iic+1]);
xoff[ii] = (fcen-ccen)/(cvc[iic+1]-cvc[iic]);
}

Expand All @@ -51,28 +51,28 @@ compute_slopes (const Dim3& lo, const Dim3& hi,
{
AMREX_PRAGMA_SIMD
for (int i = lo.x; i <= hi.x; ++i) {
slopes(i,0,0,ns) = 0.5*(u(i+1,0,0,nu)-u(i-1,0,0,nu));
slopes(i,0,0,ns) = 0.5_rt*(u(i+1,0,0,nu)-u(i-1,0,0,nu));
}

if (lo.x == slo.x && (bc.lo(0) == BCType::ext_dir || bc.lo(0) == BCType::hoextrap))
{
const int i = slo.x;
if (shi.x-slo.x >= 1) {
slopes(i,0,0,ns) = -(16./15.)*u(i-1,0,0,nu) + 0.5*u(i,0,0,nu)
+ (2./3.)*u(i+1,0,0,nu) - 0.1*u(i+2,0,0,nu);
slopes(i,0,0,ns) = -(16._rt/15._rt)*u(i-1,0,0,nu) + 0.5_rt*u(i,0,0,nu)
+ (2._rt/3._rt)*u(i+1,0,0,nu) - 0.1_rt*u(i+2,0,0,nu);
} else {
slopes(i,0,0,ns) = 0.25*(u(i+1,0,0,nu)+5.*u(i,0,0,nu)-6.*u(i-1,0,0,nu));
slopes(i,0,0,ns) = 0.25_rt*(u(i+1,0,0,nu)+5._rt*u(i,0,0,nu)-6._rt*u(i-1,0,0,nu));
}
}

if (hi.x == shi.x && (bc.hi(0) == BCType::ext_dir || bc.hi(0) == BCType::hoextrap))
{
const int i = shi.x;
if (shi.x-slo.x >= 1) {
slopes(i,0,0,ns) = (16./15.)*u(i+1,0,0,nu) - 0.5*u(i,0,0,nu)
- (2./3.)*u(i-1,0,0,nu) + 0.1*u(i-2,0,0,nu);
slopes(i,0,0,ns) = (16._rt/15._rt)*u(i+1,0,0,nu) - 0.5_rt*u(i,0,0,nu)
- (2._rt/3._rt)*u(i-1,0,0,nu) + 0.1_rt*u(i-2,0,0,nu);
} else {
slopes(i,0,0,ns) = -0.25*(u(i-1,0,0,nu)+5.*u(i,0,0,nu)-6.*u(i+1,0,0,nu));
slopes(i,0,0,ns) = -0.25_rt*(u(i-1,0,0,nu)+5._rt*u(i,0,0,nu)-6._rt*u(i+1,0,0,nu));
}
}
}
Expand All @@ -94,7 +94,7 @@ cellconslin_slopes_linlim (Box const& bx, Array4<Real> const& slopes,

AMREX_PRAGMA_SIMD
for (int i = lo.x; i <= hi.x; ++i) {
sf(i,0,0) = 1.0;
sf(i,0,0) = 1.0_rt;
}

for (int n = 0; n < ncomp; ++n)
Expand All @@ -105,14 +105,14 @@ cellconslin_slopes_linlim (Box const& bx, Array4<Real> const& slopes,
AMREX_PRAGMA_SIMD
for (int i = lo.x; i <= hi.x; ++i) {
Real cen = slopes(i,0,0,n);
Real forw = 2.0*(u(i+1,0,0,nu)-u(i ,0,0,nu));
Real back = 2.0*(u(i ,0,0,nu)-u(i-1,0,0,nu));
Real slp = (forw*back >= 0.0) ? amrex::min(amrex::Math::abs(forw),amrex::Math::abs(back)) : 0.0;
slopes(i,0,0,n) = amrex::Math::copysign(1.0,cen)*amrex::min(slp,amrex::Math::abs(cen));
if (cen != 0.0) {
Real forw = 2.0_rt*(u(i+1,0,0,nu)-u(i ,0,0,nu));
Real back = 2.0_rt*(u(i ,0,0,nu)-u(i-1,0,0,nu));
Real slp = (forw*back >= 0.0_rt) ? amrex::min(amrex::Math::abs(forw),amrex::Math::abs(back)) : 0.0_rt;
slopes(i,0,0,n) = amrex::Math::copysign(1.0_rt,cen)*amrex::min(slp,amrex::Math::abs(cen));
if (cen != 0.0_rt) {
sf(i,0,0) = amrex::min(sf(i,0,0), slopes(i,0,0,n)/cen);
} else {
sf(i,0,0) = 0.0;
sf(i,0,0) = 0.0_rt;
}
}
}
Expand Down Expand Up @@ -181,10 +181,10 @@ cellconslin_slopes_mclim (Box const& bx, Array4<Real> const& slopes,
AMREX_PRAGMA_SIMD
for (int i = lo.x; i <= hi.x; ++i) {
Real cen = slopes(i,0,0,n);
Real forw = 2.0*(u(i+1,0,0,nu)-u(i ,0,0,nu));
Real back = 2.0*(u(i ,0,0,nu)-u(i-1,0,0,nu));
Real slp = (forw*back >= 0.0) ? amrex::min(amrex::Math::abs(forw),amrex::Math::abs(back)) : 0.0;
slopes(i,0,0,n) = amrex::Math::copysign(1.0,cen)*amrex::min(slp,amrex::Math::abs(cen));
Real forw = 2.0_rt*(u(i+1,0,0,nu)-u(i ,0,0,nu));
Real back = 2.0_rt*(u(i ,0,0,nu)-u(i-1,0,0,nu));
Real slp = (forw*back >= 0.0_rt) ? amrex::min(amrex::Math::abs(forw),amrex::Math::abs(back)) : 0.0_rt;
slopes(i,0,0,n) = amrex::Math::copysign(1.0_rt,cen)*amrex::min(slp,amrex::Math::abs(cen));
}
}
}
Expand All @@ -208,12 +208,12 @@ cellconslin_fine_alpha (Box const& bx, Array4<Real> const& alpha,
const int ic = amrex::coarsen(i,ratio[0]);
const Real dummy_fine = xoff[i-vlo.x]*slopes(ic,0,0,n);

if (dummy_fine > mm(ic,0,0,n+ncomp) && dummy_fine != 0.0) {
if (dummy_fine > mm(ic,0,0,n+ncomp) && dummy_fine != 0.0_rt) {
alpha(i,0,0,n) = mm(ic,0,0,n+ncomp) / dummy_fine;
} else if (dummy_fine < mm(ic,0,0,n) && dummy_fine != 0.0) {
} else if (dummy_fine < mm(ic,0,0,n) && dummy_fine != 0.0_rt) {
alpha(i,0,0,n) = mm(ic,0,0,n) / dummy_fine;
} else {
alpha(i,0,0,n) = 1.0;
alpha(i,0,0,n) = 1.0_rt;
}
}
}
Expand All @@ -230,7 +230,7 @@ cellconslin_slopes_mmlim (Box const& bx, Array4<Real> const& slopes,
for (int n = 0; n < ncomp; ++n) {
for (int i = lo.x; i <= hi.x; ++i) {
const int ii = i*ratio[0];
Real a = 1.0;
Real a = 1.0_rt;
for (int ioff = 0; ioff < ratio[0]; ++ioff) {
a = amrex::min(a, alpha(ii+ioff,0,0,n));
}
Expand Down Expand Up @@ -265,7 +265,7 @@ nodebilin_slopes (Box const& bx, Array4<T> const& slope, Array4<T const> const&
const auto lo = amrex::lbound(bx);
const auto hi = amrex::ubound(bx);

const Real rx = 1.0/ratio[0];
const Real rx = 1.0_rt/ratio[0];

for (int n = 0; n < ncomp; ++n) {
AMREX_PRAGMA_SIMD
Expand Down
92 changes: 46 additions & 46 deletions Src/AmrCore/AMReX_Interp_2D_C.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,8 +33,8 @@ ccinterp_compute_voff (Box const& cbx, IntVect const& ratio, Geometry const& cge
const int ic = amrex::coarsen(i, ratio[0]);
const int ii = i - flo.x;
const int iic = ic - clo.x;
const Real fcen = 0.5*(fvc[ii ]+fvc[ii +1]);
const Real ccen = 0.5*(cvc[iic]+cvc[iic+1]);
const Real fcen = 0.5_rt*(fvc[ii ]+fvc[ii +1]);
const Real ccen = 0.5_rt*(cvc[iic]+cvc[iic+1]);
xoff[ii] = (fcen-ccen)/(cvc[iic+1]-cvc[iic]);
}

Expand All @@ -47,8 +47,8 @@ ccinterp_compute_voff (Box const& cbx, IntVect const& ratio, Geometry const& cge
const int jc = amrex::coarsen(j, ratio[1]);
const int jj = j - flo.y;
const int jjc = jc - clo.y;
const Real fcen = 0.5*(fvc[jj ]+fvc[jj +1]);
const Real ccen = 0.5*(cvc[jjc]+cvc[jjc+1]);
const Real fcen = 0.5_rt*(fvc[jj ]+fvc[jj +1]);
const Real ccen = 0.5_rt*(cvc[jjc]+cvc[jjc+1]);
yoff[jj] = (fcen-ccen)/(cvc[jjc+1]-cvc[jjc]);
}

Expand All @@ -66,8 +66,8 @@ compute_slopes (const Dim3& lo, const Dim3& hi,
for (int j = lo.y; j <= hi.y; ++j) {
AMREX_PRAGMA_SIMD
for (int i = lo.x; i <= hi.x; ++i) {
slopes(i,j,0,ns ) = 0.5*(u(i+1,j,0,nu)-u(i-1,j,0,nu));
slopes(i,j,0,ns+ncomp) = 0.5*(u(i,j+1,0,nu)-u(i,j-1,0,nu));
slopes(i,j,0,ns ) = 0.5_rt*(u(i+1,j,0,nu)-u(i-1,j,0,nu));
slopes(i,j,0,ns+ncomp) = 0.5_rt*(u(i,j+1,0,nu)-u(i,j-1,0,nu));
}
}

Expand All @@ -76,12 +76,12 @@ compute_slopes (const Dim3& lo, const Dim3& hi,
const int i = slo.x;
if (shi.x-slo.x >= 1) {
for (int j = lo.y; j <= hi.y; ++j) {
slopes(i,j,0,ns) = -(16./15.)*u(i-1,j,0,nu) + 0.5*u(i,j,0,nu)
+ (2./3.)*u(i+1,j,0,nu) - 0.1*u(i+2,j,0,nu);
slopes(i,j,0,ns) = -(16._rt/15._rt)*u(i-1,j,0,nu) + 0.5_rt*u(i,j,0,nu)
+ (2._rt/3._rt)*u(i+1,j,0,nu) - 0.1_rt*u(i+2,j,0,nu);
}
} else {
for (int j = lo.y; j <= hi.y; ++j) {
slopes(i,j,0,ns) = 0.25*(u(i+1,j,0,nu)+5.*u(i,j,0,nu)-6.*u(i-1,j,0,nu));
slopes(i,j,0,ns) = 0.25_rt*(u(i+1,j,0,nu)+5._rt*u(i,j,0,nu)-6._rt*u(i-1,j,0,nu));
}
}
}
Expand All @@ -91,12 +91,12 @@ compute_slopes (const Dim3& lo, const Dim3& hi,
const int i = shi.x;
if (shi.x-slo.x >= 1) {
for (int j = lo.y; j <= hi.y; ++j) {
slopes(i,j,0,ns) = (16./15.)*u(i+1,j,0,nu) - 0.5*u(i,j,0,nu)
- (2./3.)*u(i-1,j,0,nu) + 0.1*u(i-2,j,0,nu);
slopes(i,j,0,ns) = (16._rt/15._rt)*u(i+1,j,0,nu) - 0.5_rt*u(i,j,0,nu)
- (2._rt/3._rt)*u(i-1,j,0,nu) + 0.1_rt*u(i-2,j,0,nu);
}
} else {
for (int j = lo.y; j <= hi.y; ++j) {
slopes(i,j,0,ns) = -0.25*(u(i-1,j,0,nu)+5.*u(i,j,0,nu)-6.*u(i+1,j,0,nu));
slopes(i,j,0,ns) = -0.25_rt*(u(i-1,j,0,nu)+5._rt*u(i,j,0,nu)-6._rt*u(i+1,j,0,nu));
}
}
}
Expand All @@ -107,13 +107,13 @@ compute_slopes (const Dim3& lo, const Dim3& hi,
if (shi.y-slo.y >= 1) {
AMREX_PRAGMA_SIMD
for (int i = lo.x; i <= hi.x; ++i) {
slopes(i,j,0,ns+ncomp) = -(16./15.)*u(i,j-1,0,nu) + 0.5*u(i,j,0,nu)
+ (2./3.)*u(i,j+1,0,nu) - 0.1*u(i,j+2,0,nu);
slopes(i,j,0,ns+ncomp) = -(16._rt/15._rt)*u(i,j-1,0,nu) + 0.5_rt*u(i,j,0,nu)
+ (2._rt/3._rt)*u(i,j+1,0,nu) - 0.1_rt*u(i,j+2,0,nu);
}
} else {
AMREX_PRAGMA_SIMD
for (int i = lo.x; i <= hi.x; ++i) {
slopes(i,j,0,ns+ncomp) = 0.25*(u(i,j+1,0,nu)+5.*u(i,j,0,nu)-6.*u(i,j-1,0,nu));
slopes(i,j,0,ns+ncomp) = 0.25_rt*(u(i,j+1,0,nu)+5._rt*u(i,j,0,nu)-6._rt*u(i,j-1,0,nu));
}
}
}
Expand All @@ -124,13 +124,13 @@ compute_slopes (const Dim3& lo, const Dim3& hi,
if (shi.y-slo.y >= 1) {
AMREX_PRAGMA_SIMD
for (int i = lo.x; i <= hi.x; ++i) {
slopes(i,j,0,ns+ncomp) = (16./15.)*u(i,j+1,0,nu) - 0.5*u(i,j,0,nu)
- (2./3.)*u(i,j-1,0,nu) + 0.1*u(i,j-2,0,nu);
slopes(i,j,0,ns+ncomp) = (16._rt/15._rt)*u(i,j+1,0,nu) - 0.5_rt*u(i,j,0,nu)
- (2._rt/3._rt)*u(i,j-1,0,nu) + 0.1_rt*u(i,j-2,0,nu);
}
} else {
AMREX_PRAGMA_SIMD
for (int i = lo.x; i <= hi.x; ++i) {
slopes(i,j,0,ns+ncomp) = -0.25*(u(i,j-1,0,nu)+5.*u(i,j,0,nu)-6.*u(i,j+1,0,nu));
slopes(i,j,0,ns+ncomp) = -0.25_rt*(u(i,j-1,0,nu)+5._rt*u(i,j,0,nu)-6._rt*u(i,j+1,0,nu));
}
}
}
Expand All @@ -154,8 +154,8 @@ cellconslin_slopes_linlim (Box const& bx, Array4<Real> const& slopes,
for (int j = lo.y; j <= hi.y; ++j) {
AMREX_PRAGMA_SIMD
for (int i = lo.x; i <= hi.x; ++i) {
sf(i,j,0,0) = 1.0;
sf(i,j,0,1) = 1.0;
sf(i,j,0,0) = 1.0_rt;
sf(i,j,0,1) = 1.0_rt;
}
}

Expand All @@ -168,25 +168,25 @@ cellconslin_slopes_linlim (Box const& bx, Array4<Real> const& slopes,
AMREX_PRAGMA_SIMD
for (int i = lo.x; i <= hi.x; ++i) {
Real cen = slopes(i,j,0,n);
Real forw = 2.0*(u(i+1,j,0,nu)-u(i ,j,0,nu));
Real back = 2.0*(u(i ,j,0,nu)-u(i-1,j,0,nu));
Real slp = (forw*back >= 0.0) ? amrex::min(amrex::Math::abs(forw),amrex::Math::abs(back)) : 0.0;
slopes(i,j,0,n) = amrex::Math::copysign(1.0,cen)*amrex::min(slp,amrex::Math::abs(cen));
if (cen != 0.0) {
Real forw = 2.0_rt*(u(i+1,j,0,nu)-u(i ,j,0,nu));
Real back = 2.0_rt*(u(i ,j,0,nu)-u(i-1,j,0,nu));
Real slp = (forw*back >= 0.0_rt) ? amrex::min(amrex::Math::abs(forw),amrex::Math::abs(back)) : 0.0_rt;
slopes(i,j,0,n) = amrex::Math::copysign(1.0_rt,cen)*amrex::min(slp,amrex::Math::abs(cen));
if (cen != 0.0_rt) {
sf(i,j,0,0) = amrex::min(sf(i,j,0,0), slopes(i,j,0,n)/cen);
} else {
sf(i,j,0,0) = 0.0;
sf(i,j,0,0) = 0.0_rt;
}

cen = slopes(i,j,0,n+ncomp);
forw = 2.0*(u(i,j+1,0,nu)-u(i,j ,0,nu));
back = 2.0*(u(i,j ,0,nu)-u(i,j-1,0,nu));
slp = (forw*back >= 0.0) ? amrex::min(amrex::Math::abs(forw),amrex::Math::abs(back)) : 0.0;
slopes(i,j,0,n+ncomp) = amrex::Math::copysign(1.0,cen)*amrex::min(slp,amrex::Math::abs(cen));
if (cen != 0.0) {
forw = 2.0_rt*(u(i,j+1,0,nu)-u(i,j ,0,nu));
back = 2.0_rt*(u(i,j ,0,nu)-u(i,j-1,0,nu));
slp = (forw*back >= 0.0_rt) ? amrex::min(amrex::Math::abs(forw),amrex::Math::abs(back)) : 0.0_rt;
slopes(i,j,0,n+ncomp) = amrex::Math::copysign(1.0_rt,cen)*amrex::min(slp,amrex::Math::abs(cen));
if (cen != 0.0_rt) {
sf(i,j,0,1) = amrex::min(sf(i,j,0,1), slopes(i,j,0,n+ncomp)/cen);
} else {
sf(i,j,0,1) = 0.0;
sf(i,j,0,1) = 0.0_rt;
}
}
}
Expand Down Expand Up @@ -270,16 +270,16 @@ cellconslin_slopes_mclim (Box const& bx, Array4<Real> const& slopes,
AMREX_PRAGMA_SIMD
for (int i = lo.x; i <= hi.x; ++i) {
Real cen = slopes(i,j,0,n);
Real forw = 2.0*(u(i+1,j,0,nu)-u(i ,j,0,nu));
Real back = 2.0*(u(i ,j,0,nu)-u(i-1,j,0,nu));
Real slp = (forw*back >= 0.0) ? amrex::min(amrex::Math::abs(forw),amrex::Math::abs(back)) : 0.0;
slopes(i,j,0,n) = amrex::Math::copysign(1.0,cen)*amrex::min(slp,amrex::Math::abs(cen));
Real forw = 2.0_rt*(u(i+1,j,0,nu)-u(i ,j,0,nu));
Real back = 2.0_rt*(u(i ,j,0,nu)-u(i-1,j,0,nu));
Real slp = (forw*back >= 0.0_rt) ? amrex::min(amrex::Math::abs(forw),amrex::Math::abs(back)) : 0.0_rt;
slopes(i,j,0,n) = amrex::Math::copysign(1.0_rt,cen)*amrex::min(slp,amrex::Math::abs(cen));

cen = slopes(i,j,0,n+ncomp);
forw = 2.0*(u(i,j+1,0,nu)-u(i,j ,0,nu));
back = 2.0*(u(i,j ,0,nu)-u(i,j-1,0,nu));
slp = (forw*back >= 0.0) ? amrex::min(amrex::Math::abs(forw),amrex::Math::abs(back)) : 0.0;
slopes(i,j,0,n+ncomp) = amrex::Math::copysign(1.0,cen)*amrex::min(slp,amrex::Math::abs(cen));
forw = 2.0_rt*(u(i,j+1,0,nu)-u(i,j ,0,nu));
back = 2.0_rt*(u(i,j ,0,nu)-u(i,j-1,0,nu));
slp = (forw*back >= 0.0_rt) ? amrex::min(amrex::Math::abs(forw),amrex::Math::abs(back)) : 0.0_rt;
slopes(i,j,0,n+ncomp) = amrex::Math::copysign(1.0_rt,cen)*amrex::min(slp,amrex::Math::abs(cen));
}
}
}
Expand Down Expand Up @@ -309,12 +309,12 @@ cellconslin_fine_alpha (Box const& bx, Array4<Real> const& alpha,
const Real dummy_fine = xoff[i-vlo.x]*slopes(ic,jc,0,n)
+ yoff[j-vlo.y]*slopes(ic,jc,0,n+ncomp);

if (dummy_fine > mm(ic,jc,0,n+ncomp) && dummy_fine != 0.0) {
if (dummy_fine > mm(ic,jc,0,n+ncomp) && dummy_fine != 0.0_rt) {
alpha(i,j,0,n) = mm(ic,jc,0,n+ncomp) / dummy_fine;
} else if (dummy_fine < mm(ic,jc,0,n) && dummy_fine != 0.0) {
} else if (dummy_fine < mm(ic,jc,0,n) && dummy_fine != 0.0_rt) {
alpha(i,j,0,n) = mm(ic,jc,0,n) / dummy_fine;
} else {
alpha(i,j,0,n) = 1.0;
alpha(i,j,0,n) = 1.0_rt;
}
}
}
Expand All @@ -334,7 +334,7 @@ cellconslin_slopes_mmlim (Box const& bx, Array4<Real> const& slopes,
const int jj = j*ratio[1];
for (int i = lo.x; i <= hi.x; ++i) {
const int ii = i*ratio[0];
Real a = 1.0;
Real a = 1.0_rt;
for (int joff = 0; joff < ratio[1]; ++joff) {
for (int ioff = 0; ioff < ratio[0]; ++ioff) {
a = amrex::min(a, alpha(ii+ioff,jj+joff,0,n));
Expand Down Expand Up @@ -382,8 +382,8 @@ nodebilin_slopes (Box const& bx, Array4<T> const& slope, Array4<T const> const&
const auto lo = amrex::lbound(bx);
const auto hi = amrex::ubound(bx);

const Real rx = 1.0/ratio[0];
const Real ry = 1.0/ratio[1];
const Real rx = 1.0_rt/ratio[0];
const Real ry = 1.0_rt/ratio[1];

for (int n = 0; n < ncomp; ++n) {
for (int j = lo.y; j <= hi.y; ++j) {
Expand Down
Loading

0 comments on commit a017d59

Please sign in to comment.