Skip to content

Commit

Permalink
can choose number of interpolation points
Browse files Browse the repository at this point in the history
  • Loading branch information
@mh-guo
mh-guo committed Dec 6, 2024
1 parent cc9ec45 commit 06a8aa3
Show file tree
Hide file tree
Showing 2 changed files with 29 additions and 26 deletions.
52 changes: 27 additions & 25 deletions src/geodesic-grid/spherical_grid.cpp
View file
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
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

0 comments on commit 06a8aa3

Please sign in to comment.