Skip to content

Commit

Permalink
Fix multiple tracers
Browse files Browse the repository at this point in the history
  • Loading branch information
Sihe Chen authored and Sihe Chen committed Jul 12, 2024
1 parent cce2e7b commit b8990d1
Show file tree
Hide file tree
Showing 9 changed files with 81 additions and 40 deletions.
6 changes: 4 additions & 2 deletions cmake/athena.cmake
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,10 +11,12 @@ set(patch_command
${CMAKE_CURRENT_SOURCE_DIR}/patches/26.scalars_flux.patch
${CMAKE_CURRENT_SOURCE_DIR}/patches/26.adiabatic_hydro.patch
${CMAKE_CURRENT_SOURCE_DIR}/patches/27.coordinates.patch
${CMAKE_CURRENT_SOURCE_DIR}/patches/28.bvals_cc.patch
${CMAKE_CURRENT_SOURCE_DIR}/patches/29.mesh.patch
${CMAKE_CURRENT_SOURCE_DIR}/patches/30.bvals_base.patch
${CMAKE_CURRENT_SOURCE_DIR}/patches/31.task_list.patch)
${CMAKE_CURRENT_SOURCE_DIR}/patches/31.task_list.patch
${CMAKE_CURRENT_SOURCE_DIR}/patches/32.scalars.patch
${CMAKE_CURRENT_SOURCE_DIR}/patches/33.bvals_cc_cpp.patch
${CMAKE_CURRENT_SOURCE_DIR}/patches/34.bvals_cc_hpp.patch)

FetchContent_Declare(
athenapp
Expand Down
1 change: 1 addition & 0 deletions cmake/examples/exo2.cmake
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@ endmacro()

# athena variables
set_if_empty(NUMBER_GHOST_CELLS 3)
set_if_empty (NTRACER 3)

# canoe configure
set(CUBED_SPHERE ON)
Expand Down
2 changes: 2 additions & 0 deletions cmake/examples/exo3.cmake
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,3 +15,5 @@ set(NETCDF ON)
set(MPI ON)
set(PNETCDF ON)
set(RSOLVER hllc_transform)

set (NTRACER 3)
12 changes: 11 additions & 1 deletion examples/2023-Chen-exo3/steady_zonal.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@
#include <athena/hydro/hydro.hpp>
#include <athena/mesh/mesh.hpp>
#include <athena/parameter_input.hpp>
#include <athena/scalars/scalars.hpp>

// application
#include <application/application.hpp>
Expand Down Expand Up @@ -80,6 +81,11 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) {
pexo3->GetLatLon(&lat, &lon, k, j, i);
phydro->w(IDN, k, j, i) =
g * h0 - (a * om_earth * u0 + u0 * u0 / 2) * sin(lat) * sin(lat);
if ((lat > PI / 2.0 * 0.9) && (lat < PI / 2.0))
pscalars->s(2, k, j, i) = 1.0;
if ((lat > PI / 2.0 * 0.5) && (lat < PI / 2.0 * 0.6))
pscalars->s(1, k, j, i) = 1.0;

Real U = u0 * cos(lat);
Real V = 0.0;
Real Vy, Vz;
Expand Down Expand Up @@ -118,16 +124,20 @@ void MeshBlock::UserWorkBeforeOutput(ParameterInput *pin) {
user_out_var(2, k, j, i) = U;
user_out_var(3, k, j, i) = V;
user_out_var(4, k, j, i) = delta / (C * D);
user_out_var(5, k, j, i) = pscalars->s(1, k, j, i);
user_out_var(6, k, j, i) = pscalars->s(2, k, j, i);
}
}

void MeshBlock::InitUserMeshBlockData(ParameterInput *pin) {
AllocateUserOutputVariables(5);
AllocateUserOutputVariables(7);
SetUserOutputVariableName(0, "lat");
SetUserOutputVariableName(1, "lon");
SetUserOutputVariableName(2, "U");
SetUserOutputVariableName(3, "V");
SetUserOutputVariableName(4, "sqrtg");
SetUserOutputVariableName(5, "s1");
SetUserOutputVariableName(6, "s2");
}

void Mesh::InitUserMeshData(ParameterInput *pin) {
Expand Down
12 changes: 12 additions & 0 deletions patches/32.scalars.patch
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,12 @@
diff --git a/src/scalars/scalars.cpp b/src/scalars/scalars.cpp
index 8471fb20..2f3ad1ce 100644
--- a/src/scalars/scalars.cpp
+++ b/src/scalars/scalars.cpp
@@ -84,6 +84,7 @@ PassiveScalars::PassiveScalars(MeshBlock *pmb, ParameterInput *pin) :

// enroll CellCenteredBoundaryVariable object
sbvar.bvar_index = pmb->pbval->bvars.size();
+ sbvar.TypeFlag = 1; // PassiveScalars
pmb->pbval->bvars.push_back(&sbvar);
pmb->pbval->bvars_main_int.push_back(&sbvar);
if (STS_ENABLED) {
14 changes: 7 additions & 7 deletions patches/28.bvals_cc.patch → patches/33.bvals_cc_cpp.patch
View file
Original file line number Diff line number Diff line change
@@ -1,24 +1,24 @@
diff --git a/src/bvals/cc/bvals_cc.cpp b/src/bvals/cc/bvals_cc.cpp
index eb6f6039..d8b328b0 100644
index 223962de..74b3ab79 100644
--- a/src/bvals/cc/bvals_cc.cpp
+++ b/src/bvals/cc/bvals_cc.cpp
@@ -30,6 +30,9 @@
#include "../../utils/buffer_utils.hpp"
@@ -31,6 +31,9 @@
#include "../bvals.hpp"
#include "bvals_cc.hpp"
+
+#include <configure.hpp>
+#include <exo3/cubed_sphere_utility.hpp>
+
// MPI header
#ifdef MPI_PARALLEL
@@ -297,7 +298,12 @@ int CellCenteredBoundaryVariable::LoadBoundaryBufferSameLevel(Real *buf,
#include <mpi.h>
@@ -297,7 +300,12 @@ int CellCenteredBoundaryVariable::LoadBoundaryBufferSameLevel(Real *buf,
ek = (nb.ni.ox3 < 0) ? (pmb->ks + NGHOST - 1) : pmb->ke;
int p = 0;
AthenaArray<Real> &var = *var_cc;
+#ifdef CUBED_SPHERE
+// nl_, nu_, after buf var
+ CubedSphereUtility::PackData(var, buf, nl_, nu_, si, ei, sj, ej, sk, ek, p, nb.ni.ox1, nb.ni.ox2, nb.ni.ox3, pmb->loc);
+ CubedSphereUtility::PackData(var, buf, nl_, nu_, si, ei, sj, ej, sk, ek, p, nb.ni.ox1, nb.ni.ox2, nb.ni.ox3, pmb->loc, TypeFlag);
+#else
BufferUtility::PackData(var, buf, nl_, nu_, si, ei, sj, ej, sk, ek, p);
+#endif
Expand Down
14 changes: 14 additions & 0 deletions patches/34.bvals_cc_hpp.patch
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,14 @@
diff --git a/src/bvals/cc/bvals_cc.hpp b/src/bvals/cc/bvals_cc.hpp
index d167d729..ce2f8dcc 100644
--- a/src/bvals/cc/bvals_cc.hpp
+++ b/src/bvals/cc/bvals_cc.hpp
@@ -39,6 +39,9 @@ class CellCenteredBoundaryVariable : public BoundaryVariable {
AthenaArray<Real> *var_flux, bool fflux, int flag);
~CellCenteredBoundaryVariable();

+ // Type information
+ int TypeFlag = 2; // Hydro by default
+
//! \note
//! may want to rebind var_cc to u,u1,u2,w,w1, etc. registers for time integrator logic.
//! Also, derived class HydroBoundaryVariable needs to keep switching var and coarse_var
58 changes: 29 additions & 29 deletions src/exo3/cubed_sphere_utility.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -99,7 +99,7 @@ Real CalculateInterpLocations(int loc_n, int N_blk, int k, bool GhostZone) {
void InteprolateX2(const AthenaArray<Real> &src, AthenaArray<Real> &tgt,
LogicalLocation const &loc, int DirInv, int TgtSide,
int TgtID, int sn, int en, int si, int ei, int sj, int ej,
int sk, int ek) {
int sk, int ek, int TypeFlag) {
// Interpolation along X2 (j) axis, used before sending data to X3 (k) axis
// Get the local indices
int lv2_lx2 = loc.lx2 >> (loc.level - 2);
Expand Down Expand Up @@ -159,7 +159,7 @@ void InteprolateX2(const AthenaArray<Real> &src, AthenaArray<Real> &tgt,
Real y1 = src_x2[src_pointer];
Real y2 = src_x2[src_pointer + 1];
Real yq = tgt_x2[j - sj];
if (n == IVY || n == IVZ) {
if ((n == IVY || n == IVZ) && (TypeFlag == 2)) {
// Projection needed, find the tgt locations first
Real v1y =
src(IVY, k, src_pointer - n_start - N_blk / 2 + sj,
Expand Down Expand Up @@ -230,7 +230,7 @@ void InteprolateX2(const AthenaArray<Real> &src, AthenaArray<Real> &tgt,
void InteprolateX3(const AthenaArray<Real> &src, AthenaArray<Real> &tgt,
LogicalLocation const &loc, int DirInv, int TgtSide,
int TgtID, int sn, int en, int si, int ei, int sj, int ej,
int sk, int ek) {
int sk, int ek, int TypeFlag) {
// Interpolation along X3 (k) axis, used before sending data to ghost zone in
// X2 (j) direction Get the local indices
int lv2_lx2 = loc.lx2 >> (loc.level - 2);
Expand Down Expand Up @@ -288,7 +288,7 @@ void InteprolateX3(const AthenaArray<Real> &src, AthenaArray<Real> &tgt,
Real y1 = src_x3[src_pointer];
Real y2 = src_x3[src_pointer + 1];
Real yq = tgt_x3[k - sk];
if (n == IVY || n == IVZ) {
if ((n == IVY || n == IVZ) && (TypeFlag == 2)){
// Projection needed, find the tgt locations first
Real v1y = src(IVY, src_pointer - n_start - N_blk / 2 + sk, j, i);
Real v1z = src(IVZ, src_pointer - n_start - N_blk / 2 + sk, j, i);
Expand Down Expand Up @@ -377,7 +377,7 @@ void InteprolateX3(const AthenaArray<Real> &src, AthenaArray<Real> &tgt,

void PackData(const AthenaArray<Real> &src, Real *buf, int sn, int en, int si,
int ei, int sj, int ej, int sk, int ek, int &offset, int ox1,
int ox2, int ox3, LogicalLocation const &loc) {
int ox2, int ox3, LogicalLocation const &loc, int TypeFlag) {
// Find the block ID
int blockID = CubedSphere::FindBlockID(loc);

Expand Down Expand Up @@ -419,31 +419,31 @@ void PackData(const AthenaArray<Real> &src, Real *buf, int sn, int en, int si,
if (ox3 == 1 && local_lx3 == bound_lim) {
InteprolateX2(src, interpolatedSrc, loc, dinv[blockID - 1][3],
tgside[blockID - 1][3], tgbid[blockID - 1][3], sn, en, si,
ei, sj, ej, sk, ek);
ei, sj, ej, sk, ek, TypeFlag);
PackDataR1(interpolatedSrc, buf, 0, en - sn, 0, ei - si, 0, ej - sj, 0,
ek - sk, offset);
return;
}
if (ox2 == 1 && local_lx2 == bound_lim) {
InteprolateX3(src, interpolatedSrc, loc, dinv[blockID - 1][1],
tgside[blockID - 1][1], tgbid[blockID - 1][1], sn, en, si,
ei, sj, ej, sk, ek);
ei, sj, ej, sk, ek, TypeFlag);
PackDataR0(interpolatedSrc, buf, 0, en - sn, 0, ei - si, 0, ej - sj, 0,
ek - sk, offset);
return;
}
if (ox2 == -1 && local_lx2 == 0) {
InteprolateX3(src, interpolatedSrc, loc, dinv[blockID - 1][0],
tgside[blockID - 1][0], tgbid[blockID - 1][0], sn, en, si,
ei, sj, ej, sk, ek);
ei, sj, ej, sk, ek, TypeFlag);
PackDataR2(interpolatedSrc, buf, 0, en - sn, 0, ei - si, 0, ej - sj, 0,
ek - sk, offset);
return;
}
if (ox3 == -1 && local_lx3 == 0) {
InteprolateX2(src, interpolatedSrc, loc, dinv[blockID - 1][2],
tgside[blockID - 1][2], tgbid[blockID - 1][2], sn, en, si,
ei, sj, ej, sk, ek);
ei, sj, ej, sk, ek, TypeFlag);
PackDataR3(interpolatedSrc, buf, 0, en - sn, 0, ei - si, 0, ej - sj, 0,
ek - sk, offset);
return;
Expand All @@ -453,31 +453,31 @@ void PackData(const AthenaArray<Real> &src, Real *buf, int sn, int en, int si,
if (ox3 == 1 && local_lx3 == bound_lim) {
InteprolateX2(src, interpolatedSrc, loc, dinv[blockID - 1][3],
tgside[blockID - 1][3], tgbid[blockID - 1][3], sn, en, si,
ei, sj, ej, sk, ek);
ei, sj, ej, sk, ek, TypeFlag);
PackDataR0(interpolatedSrc, buf, 0, en - sn, 0, ei - si, 0, ej - sj, 0,
ek - sk, offset);
return;
}
if (ox2 == 1 && local_lx2 == bound_lim) {
InteprolateX3(src, interpolatedSrc, loc, dinv[blockID - 1][1],
tgside[blockID - 1][1], tgbid[blockID - 1][1], sn, en, si,
ei, sj, ej, sk, ek);
ei, sj, ej, sk, ek, TypeFlag);
PackDataR0(interpolatedSrc, buf, 0, en - sn, 0, ei - si, 0, ej - sj, 0,
ek - sk, offset);
return;
}
if (ox2 == -1 && local_lx2 == 0) {
InteprolateX3(src, interpolatedSrc, loc, dinv[blockID - 1][0],
tgside[blockID - 1][0], tgbid[blockID - 1][0], sn, en, si,
ei, sj, ej, sk, ek);
ei, sj, ej, sk, ek, TypeFlag);
PackDataR0(interpolatedSrc, buf, 0, en - sn, 0, ei - si, 0, ej - sj, 0,
ek - sk, offset);
return;
}
if (ox3 == -1 && local_lx3 == 0) {
InteprolateX2(src, interpolatedSrc, loc, dinv[blockID - 1][2],
tgside[blockID - 1][2], tgbid[blockID - 1][2], sn, en, si,
ei, sj, ej, sk, ek);
ei, sj, ej, sk, ek, TypeFlag);
PackDataR0(interpolatedSrc, buf, 0, en - sn, 0, ei - si, 0, ej - sj, 0,
ek - sk, offset);
return;
Expand All @@ -487,31 +487,31 @@ void PackData(const AthenaArray<Real> &src, Real *buf, int sn, int en, int si,
if (ox2 == 1 && local_lx2 == bound_lim) {
InteprolateX3(src, interpolatedSrc, loc, dinv[blockID - 1][1],
tgside[blockID - 1][1], tgbid[blockID - 1][1], sn, en, si,
ei, sj, ej, sk, ek);
ei, sj, ej, sk, ek, TypeFlag);
PackDataR3(interpolatedSrc, buf, 0, en - sn, 0, ei - si, 0, ej - sj, 0,
ek - sk, offset);
return;
}
if (ox2 == -1 && local_lx2 == 0) {
InteprolateX3(src, interpolatedSrc, loc, dinv[blockID - 1][0],
tgside[blockID - 1][0], tgbid[blockID - 1][0], sn, en, si,
ei, sj, ej, sk, ek);
ei, sj, ej, sk, ek, TypeFlag);
PackDataR1(interpolatedSrc, buf, 0, en - sn, 0, ei - si, 0, ej - sj, 0,
ek - sk, offset);
return;
}
if (ox3 == 1 && local_lx3 == bound_lim) {
InteprolateX2(src, interpolatedSrc, loc, dinv[blockID - 1][3],
tgside[blockID - 1][3], tgbid[blockID - 1][3], sn, en, si,
ei, sj, ej, sk, ek);
ei, sj, ej, sk, ek, TypeFlag);
PackDataR0(interpolatedSrc, buf, 0, en - sn, 0, ei - si, 0, ej - sj, 0,
ek - sk, offset);
return;
}
if (ox3 == -1 && local_lx3 == 0) {
InteprolateX2(src, interpolatedSrc, loc, dinv[blockID - 1][2],
tgside[blockID - 1][2], tgbid[blockID - 1][2], sn, en, si,
ei, sj, ej, sk, ek);
ei, sj, ej, sk, ek, TypeFlag);
PackDataR0(interpolatedSrc, buf, 0, en - sn, 0, ei - si, 0, ej - sj, 0,
ek - sk, offset);
return;
Expand All @@ -521,31 +521,31 @@ void PackData(const AthenaArray<Real> &src, Real *buf, int sn, int en, int si,
if (ox2 == 1 && local_lx2 == bound_lim) {
InteprolateX3(src, interpolatedSrc, loc, dinv[blockID - 1][1],
tgside[blockID - 1][1], tgbid[blockID - 1][1], sn, en, si,
ei, sj, ej, sk, ek);
ei, sj, ej, sk, ek, TypeFlag);
PackDataR1(interpolatedSrc, buf, 0, en - sn, 0, ei - si, 0, ej - sj, 0,
ek - sk, offset);
return;
}
if (ox2 == -1 && local_lx2 == 0) {
InteprolateX3(src, interpolatedSrc, loc, dinv[blockID - 1][0],
tgside[blockID - 1][0], tgbid[blockID - 1][0], sn, en, si,
ei, sj, ej, sk, ek);
ei, sj, ej, sk, ek, TypeFlag);
PackDataR3(interpolatedSrc, buf, 0, en - sn, 0, ei - si, 0, ej - sj, 0,
ek - sk, offset);
return;
}
if (ox3 == 1 && local_lx3 == bound_lim) {
InteprolateX2(src, interpolatedSrc, loc, dinv[blockID - 1][3],
tgside[blockID - 1][3], tgbid[blockID - 1][3], sn, en, si,
ei, sj, ej, sk, ek);
ei, sj, ej, sk, ek, TypeFlag);
PackDataR0(interpolatedSrc, buf, 0, en - sn, 0, ei - si, 0, ej - sj, 0,
ek - sk, offset);
return;
}
if (ox3 == -1 && local_lx3 == 0) {
InteprolateX2(src, interpolatedSrc, loc, dinv[blockID - 1][2],
tgside[blockID - 1][2], tgbid[blockID - 1][2], sn, en, si,
ei, sj, ej, sk, ek);
ei, sj, ej, sk, ek, TypeFlag);
PackDataR0(interpolatedSrc, buf, 0, en - sn, 0, ei - si, 0, ej - sj, 0,
ek - sk, offset);
return;
Expand All @@ -555,31 +555,31 @@ void PackData(const AthenaArray<Real> &src, Real *buf, int sn, int en, int si,
if (ox3 == 1 && local_lx3 == bound_lim) {
InteprolateX2(src, interpolatedSrc, loc, dinv[blockID - 1][3],
tgside[blockID - 1][3], tgbid[blockID - 1][3], sn, en, si,
ei, sj, ej, sk, ek);
ei, sj, ej, sk, ek, TypeFlag);
PackDataR3(interpolatedSrc, buf, 0, en - sn, 0, ei - si, 0, ej - sj, 0,
ek - sk, offset);
return;
}
if (ox2 == 1 && local_lx2 == bound_lim) {
InteprolateX3(src, interpolatedSrc, loc, dinv[blockID - 1][1],
tgside[blockID - 1][1], tgbid[blockID - 1][1], sn, en, si,
ei, sj, ej, sk, ek);
ei, sj, ej, sk, ek, TypeFlag);
PackDataR2(interpolatedSrc, buf, 0, en - sn, 0, ei - si, 0, ej - sj, 0,
ek - sk, offset);
return;
}
if (ox3 == -1 && local_lx3 == 0) {
InteprolateX2(src, interpolatedSrc, loc, dinv[blockID - 1][2],
tgside[blockID - 1][2], tgbid[blockID - 1][2], sn, en, si,
ei, sj, ej, sk, ek);
ei, sj, ej, sk, ek, TypeFlag);
PackDataR1(interpolatedSrc, buf, 0, en - sn, 0, ei - si, 0, ej - sj, 0,
ek - sk, offset);
return;
}
if (ox2 == -1 && local_lx2 == 0) {
InteprolateX3(src, interpolatedSrc, loc, dinv[blockID - 1][0],
tgside[blockID - 1][0], tgbid[blockID - 1][0], sn, en, si,
ei, sj, ej, sk, ek);
ei, sj, ej, sk, ek, TypeFlag);
PackDataR0(interpolatedSrc, buf, 0, en - sn, 0, ei - si, 0, ej - sj, 0,
ek - sk, offset);
return;
Expand All @@ -589,31 +589,31 @@ void PackData(const AthenaArray<Real> &src, Real *buf, int sn, int en, int si,
if (ox2 == 1 && local_lx2 == bound_lim) {
InteprolateX3(src, interpolatedSrc, loc, dinv[blockID - 1][1],
tgside[blockID - 1][1], tgbid[blockID - 1][1], sn, en, si,
ei, sj, ej, sk, ek);
ei, sj, ej, sk, ek, TypeFlag);
PackDataR2(interpolatedSrc, buf, 0, en - sn, 0, ei - si, 0, ej - sj, 0,
ek - sk, offset);
return;
}
if (ox2 == -1 && local_lx2 == 0) {
InteprolateX3(src, interpolatedSrc, loc, dinv[blockID - 1][0],
tgside[blockID - 1][0], tgbid[blockID - 1][0], sn, en, si,
ei, sj, ej, sk, ek);
ei, sj, ej, sk, ek, TypeFlag);
PackDataR2(interpolatedSrc, buf, 0, en - sn, 0, ei - si, 0, ej - sj, 0,
ek - sk, offset);
return;
}
if (ox3 == 1 && local_lx3 == bound_lim) {
InteprolateX2(src, interpolatedSrc, loc, dinv[blockID - 1][3],
tgside[blockID - 1][3], tgbid[blockID - 1][3], sn, en, si,
ei, sj, ej, sk, ek);
ei, sj, ej, sk, ek, TypeFlag);
PackDataR0(interpolatedSrc, buf, 0, en - sn, 0, ei - si, 0, ej - sj, 0,
ek - sk, offset);
return;
}
if (ox3 == -1 && local_lx3 == 0) {
InteprolateX2(src, interpolatedSrc, loc, dinv[blockID - 1][2],
tgside[blockID - 1][2], tgbid[blockID - 1][2], sn, en, si,
ei, sj, ej, sk, ek);
ei, sj, ej, sk, ek, TypeFlag);
PackDataR0(interpolatedSrc, buf, 0, en - sn, 0, ei - si, 0, ej - sj, 0,
ek - sk, offset);
return;
Expand Down
Loading

0 comments on commit b8990d1

Please sign in to comment.