-
Notifications
You must be signed in to change notification settings - Fork 341
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
add geometric terms for spherical 2D support. #4141
Changes from 1 commit
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -164,6 +164,7 @@ CoordSys::UpperIndex(const Real* point) const noexcept | |
IntVect ix; | ||
for (int k = 0; k < AMREX_SPACEDIM; k++) | ||
{ | ||
// +1 for UpperIndex? | ||
ix[k] = (int) ((point[k]-offset[k])/dx[k]); | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Maybe we should do +1 to the int, not the real number before the cast. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. I moved the +1 outside int() |
||
} | ||
return ix; | ||
|
@@ -330,6 +331,8 @@ CoordSys::GetEdgeVolCoord (Vector<Real>& vc, | |
GetEdgeLoc(vc,region,dir); | ||
// | ||
// In R direction of RZ, vol coord = (r^2)/2 | ||
// In R direction of SPHERICAL, vol coord = (r^3)/3 | ||
// In theta direction of SPHERICAL, vol coord = -cos(theta) | ||
// | ||
#if (AMREX_SPACEDIM == 2) | ||
if (dir == 0 && c_sys == RZ) | ||
|
@@ -342,6 +345,29 @@ CoordSys::GetEdgeVolCoord (Vector<Real>& vc, | |
vc[i] = 0.5_rt*r*r; | ||
} | ||
} | ||
else if (c_sys == SPHERICAL) | ||
{ | ||
if (dir == 0) | ||
{ | ||
int len = static_cast<int>(vc.size()); | ||
AMREX_PRAGMA_SIMD | ||
for (int i = 0; i < len; i++) | ||
{ | ||
Real r = vc[i]; | ||
vc[i] = r*r*r/3.0_rt; | ||
} | ||
} | ||
else | ||
{ | ||
int len = static_cast<int>(vc.size()); | ||
AMREX_PRAGMA_SIMD | ||
for (int i = 0; i < len; i++) | ||
{ | ||
Real theta = vc[i]; | ||
vc[i] = -std::cos(theta); | ||
} | ||
} | ||
} | ||
#elif (AMREX_SPACEDIM == 1) | ||
if (c_sys == SPHERICAL) | ||
{ | ||
|
@@ -365,8 +391,11 @@ CoordSys::GetCellVolCoord (Vector<Real>& vc, | |
// are identical to physical distance from axis. | ||
// | ||
GetCellLoc(vc,region,dir); | ||
|
||
// | ||
// In R direction of RZ, vol coord = (r^2)/2. | ||
// In R direction of RZ, vol coord = (r^2)/2 | ||
// In R direction of SPHERICAL, vol coord = (r^3)/3 | ||
// In theta direction of SPHERICAL, vol coord = -cos(theta) | ||
// | ||
#if (AMREX_SPACEDIM == 2) | ||
if (dir == 0 && c_sys == RZ) | ||
|
@@ -379,6 +408,29 @@ CoordSys::GetCellVolCoord (Vector<Real>& vc, | |
vc[i] = 0.5_rt*r*r; | ||
} | ||
} | ||
else if (c_sys == SPHERICAL) | ||
{ | ||
if (dir == 0) | ||
{ | ||
int len = static_cast<int>(vc.size()); | ||
AMREX_PRAGMA_SIMD | ||
for (int i = 0; i < len; i++) | ||
{ | ||
Real r = vc[i]; | ||
vc[i] = r*r*r/3.0_rt; | ||
} | ||
} | ||
else | ||
{ | ||
int len = static_cast<int>(vc.size()); | ||
AMREX_PRAGMA_SIMD | ||
for (int i = 0; i < len; i++) | ||
{ | ||
Real theta = vc[i]; | ||
vc[i] = -std::cos(theta); | ||
} | ||
} | ||
} | ||
#elif (AMREX_SPACEDIM == 1) | ||
if (c_sys == SPHERICAL) { | ||
int len = static_cast<int>(vc.size()); | ||
|
@@ -461,7 +513,10 @@ CoordSys::Volume (const Real xlo[AMREX_SPACEDIM], | |
*(xhi[2]-xlo[2])); | ||
#if (AMREX_SPACEDIM==2) | ||
case RZ: | ||
return static_cast<Real>(0.5*TWOPI)*(xhi[1]-xlo[1])*(xhi[0]*xhi[0]-xlo[0]*xlo[0]); | ||
return static_cast<Real>(0.5*TWOPI)*(xhi[1]-xlo[1])*(xhi[0]+xlo[0])*(xhi[0]-xlo[0]); | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. this will cause roundoff diffs when merged There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. should I change it back? There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Yes, that would be better. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. okay, I changed it back. |
||
case SPHERICAL: | ||
return static_cast<Real>(TWOPI/3.)*(std::cos(xlo[1])-std::cos(xhi[1])) * | ||
(xhi[0]-xlo[0])*(xhi[0]*xhi[0]+xhi[0]*xlo[0]+xlo[0]*xlo[0]); | ||
#endif | ||
default: | ||
AMREX_ASSERT(0); | ||
|
@@ -496,6 +551,16 @@ CoordSys::AreaLo (const IntVect& point, int dir) const noexcept // NOLINT(readab | |
AMREX_ASSERT(0); | ||
} | ||
return 0._rt; // to silent compiler warning | ||
case SPHERICAL: | ||
LoNode(point,xlo); | ||
switch (dir) | ||
{ | ||
case 0: return Real(TWOPI)*xlo[0]*xlo[0]*(std::cos(xlo[1]) - std::cos(xlo[1]+dx[1])); | ||
case 1: return (xlo[0]+xlo[0]+dx[0])*dx[0]*std::sin(xlo[1])*static_cast<Real>(0.5*TWOPI); | ||
default: | ||
AMREX_ASSERT(0); | ||
} | ||
return 0._rt; // to silent compiler warning | ||
default: | ||
AMREX_ASSERT(0); | ||
} | ||
|
@@ -540,6 +605,16 @@ CoordSys::AreaHi (const IntVect& point, int dir) const noexcept // NOLINT(readab | |
AMREX_ASSERT(0); | ||
} | ||
return 0._rt; // to silent compiler warning | ||
case SPHERICAL: | ||
HiNode(point,xhi); | ||
switch (dir) | ||
{ | ||
case 0: return Real(TWOPI)*xhi[0]*xhi[0]*(std::cos(xhi[1]-dx[1]) - std::cos(xhi[1])); | ||
case 1: return (xhi[0]+xhi[0]-dx[0])*dx[0]*std::sin(xhi[1])*static_cast<Real>(0.5*TWOPI); | ||
default: | ||
AMREX_ASSERT(0); | ||
} | ||
return 0._rt; // to silent compiler warning | ||
default: | ||
AMREX_ASSERT(0); | ||
} | ||
|
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
btw, I was wondering if
UpperIndex
should have a +1There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
it looks like it should, but it also appears that this isn't used anywhere