Skip to content
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.

Sign up for GitHub

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

Now can choose number of interpolation points #621

Draft
wants to merge 1 commit into
base: main
Choose a base branch
from
Draft

Now can choose number of interpolation points #621

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
52 changes: 27 additions & 25 deletions src/geodesic-grid/spherical_grid.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@
//----------------------------------------------------------------------------------------
// constructor, initializes data structures and parameters

SphericalGrid::SphericalGrid(MeshBlockPack *ppack, int nlev, Real rad):
SphericalGrid::SphericalGrid(MeshBlockPack *ppack, int nlev, Real rad, int nintp):
GeodesicGrid(nlev,true,false),
pmy_pack(ppack),
radius(rad),
Expand All @@ -32,10 +32,10 @@ SphericalGrid::SphericalGrid(MeshBlockPack *ppack, int nlev, Real rad):
interp_wghts("interp_wghts",1,1,1),
interp_vals("interp_vals",1,1) {
// reallocate and set interpolation coordinates, indices, and weights
int &ng = pmy_pack->pmesh->mb_indcs.ng;
ninterp = (nintp <= 0) ? pmy_pack->pmesh->mb_indcs.ng*2 : nintp;
Kokkos::realloc(interp_coord,nangles,3);
Kokkos::realloc(interp_indcs,nangles,4);
Kokkos::realloc(interp_wghts,nangles,2*ng,3);
Kokkos::realloc(interp_wghts,nangles,ninterp,3);

// Call functions to prepare SphericalGrid object for interpolation
SetInterpolationCoordinates();
Expand Down Expand Up @@ -120,17 +120,18 @@ void SphericalGrid::SetInterpolationIndices() {
Real &dx3 = size.h_view(m).dx3;

// save MeshBlock and zone indicies for nearest position to spherical patch center
// if this angle position resides in this MeshBlock
if ((rcoord.h_view(n,0) >= x1min && rcoord.h_view(n,0) <= x1max) &&
(rcoord.h_view(n,1) >= x2min && rcoord.h_view(n,1) <= x2max) &&
(rcoord.h_view(n,2) >= x3min && rcoord.h_view(n,2) <= x3max)) {
// if this angle position resides in the active cells of this MeshBlock without
// double count between different ranks
if ((rcoord.h_view(n,0) >= x1min && rcoord.h_view(n,0) < x1max) &&
(rcoord.h_view(n,1) >= x2min && rcoord.h_view(n,1) < x2max) &&
(rcoord.h_view(n,2) >= x3min && rcoord.h_view(n,2) < x3max)) {
iindcs.h_view(n,0) = m;
iindcs.h_view(n,1) = static_cast<int>(std::floor((rcoord.h_view(n,0)-
(x1min+dx1/2.0))/dx1));
(x1min))/dx1));
iindcs.h_view(n,2) = static_cast<int>(std::floor((rcoord.h_view(n,1)-
(x2min+dx2/2.0))/dx2));
(x2min))/dx2));
iindcs.h_view(n,3) = static_cast<int>(std::floor((rcoord.h_view(n,2)-
(x3min+dx3/2.0))/dx3));
(x3min))/dx3));
}
}
}
Expand All @@ -149,7 +150,6 @@ void SphericalGrid::SetInterpolationIndices() {
void SphericalGrid::SetInterpolationWeights() {
auto &indcs = pmy_pack->pmesh->mb_indcs;
auto &size = pmy_pack->pmb->mb_size;
int &ng = indcs.ng;

auto &iindcs = interp_indcs;
auto &iwghts = interp_wghts;
Expand All @@ -161,7 +161,7 @@ void SphericalGrid::SetInterpolationWeights() {
int &ii3 = iindcs.h_view(n,3);

if (ii0==-1) { // angle not on this rank
for (int i=0; i<2*ng; ++i) {
for (int i=0; i<ninterp; ++i) {
iwghts.h_view(n,i,0) = 0.0;
iwghts.h_view(n,i,1) = 0.0;
iwghts.h_view(n,i,2) = 0.0;
Expand All @@ -181,20 +181,21 @@ void SphericalGrid::SetInterpolationWeights() {
Real &x3max = size.h_view(ii0).x3max;

// set interpolation weights
for (int i=0; i<2*ng; ++i) {
int nleft = ninterp/2;
for (int i=0; i<ninterp; ++i) {
iwghts.h_view(n,i,0) = 1.;
iwghts.h_view(n,i,1) = 1.;
iwghts.h_view(n,i,2) = 1.;
for (int j=0; j<2*ng; ++j) {
Real x1vpi1 = CellCenterX(ii1-nleft+i, indcs.nx1, x1min, x1max);
Real x2vpi1 = CellCenterX(ii2-nleft+i, indcs.nx2, x2min, x2max);
Real x3vpi1 = CellCenterX(ii3-nleft+i, indcs.nx3, x3min, x3max);
for (int j=0; j<ninterp; ++j) {
if (j != i) {
Real x1vpi1 = CellCenterX(ii1-ng+i+1, indcs.nx1, x1min, x1max);
Real x1vpj1 = CellCenterX(ii1-ng+j+1, indcs.nx1, x1min, x1max);
Real x1vpj1 = CellCenterX(ii1-nleft+j, indcs.nx1, x1min, x1max);
iwghts.h_view(n,i,0) *= (x0-x1vpj1)/(x1vpi1-x1vpj1);
Real x2vpi1 = CellCenterX(ii2-ng+i+1, indcs.nx2, x2min, x2max);
Real x2vpj1 = CellCenterX(ii2-ng+j+1, indcs.nx2, x2min, x2max);
Real x2vpj1 = CellCenterX(ii2-nleft+j, indcs.nx2, x2min, x2max);
iwghts.h_view(n,i,1) *= (y0-x2vpj1)/(x2vpi1-x2vpj1);
Real x3vpi1 = CellCenterX(ii3-ng+i+1, indcs.nx3, x3min, x3max);
Real x3vpj1 = CellCenterX(ii3-ng+j+1, indcs.nx3, x3min, x3max);
Real x3vpj1 = CellCenterX(ii3-nleft+j, indcs.nx3, x3min, x3max);
iwghts.h_view(n,i,2) *= (z0-x3vpj1)/(x3vpi1-x3vpj1);
}
}
Expand Down Expand Up @@ -223,9 +224,10 @@ void SphericalGrid::InterpolateToSphere(int nvars, DvceArray5D<Real> &val) {
// capturing variables for kernel
auto &indcs = pmy_pack->pmesh->mb_indcs;
int &is = indcs.is; int &js = indcs.js; int &ks = indcs.ks;
int &ng = indcs.ng;
int nang1 = nangles - 1;
int nvar1 = nvars - 1;
int nintp = ninterp;
int nleft = nintp/2;

// reallocate container
Kokkos::realloc(interp_vals,nangles,nvars);
Expand All @@ -244,11 +246,11 @@ void SphericalGrid::InterpolateToSphere(int nvars, DvceArray5D<Real> &val) {
ivals.d_view(n,v) = 0.0;
} else {
Real int_value = 0.0;
for (int i=0; i<2*ng; i++) {
for (int j=0; j<2*ng; j++) {
for (int k=0; k<2*ng; k++) {
for (int i=0; i<nintp; i++) {
for (int j=0; j<nintp; j++) {
for (int k=0; k<nintp; k++) {
Real iwght = iwghts.d_view(n,i,0)*iwghts.d_view(n,j,1)*iwghts.d_view(n,k,2);
int_value += iwght*val(ii0,v,ii3-(ng-k-ks)+1,ii2-(ng-j-js)+1,ii1-(ng-i-is)+1);
int_value += iwght*val(ii0,v,ii3+k+ks-nleft,ii2+j+js-nleft,ii1+i+is-nleft);
}
}
}
Expand Down
3 changes: 2 additions & 1 deletion src/geodesic-grid/spherical_grid.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,10 +21,11 @@ class MeshBlockPack;
class SphericalGrid: public GeodesicGrid {
public:
// Creates a geodesic grid with refinement level nlev and radius rad
SphericalGrid(MeshBlockPack *pmy_pack, int nlev, Real rad);
SphericalGrid(MeshBlockPack *pmy_pack, int nlev, Real rad, int ninterp = -1);
~SphericalGrid();

Real radius; // const radius for SphericalGrid
int ninterp; // number of interpolation points along each dimension
DualArray2D<Real> interp_coord; // Cartesian coordinates for grid points
DualArray2D<Real> interp_vals; // container for data interpolated to sphere
void InterpolateToSphere(int nvars, DvceArray5D<Real> &val); // interpolate to sphere
Expand Down