From 7cc524b78b2807f520dc4ad2c4f92329a1a83209 Mon Sep 17 00:00:00 2001 From: Miikka Vaisala Date: Wed, 7 Aug 2019 14:57:51 +0800 Subject: [PATCH] Adapting for autotest but i, j, k indexing is confusing. --- acc/mhd_solver/stencil_defines.h | 4 +- src/standalone/model/model_rk3.cc | 118 +++++++++++++++--------------- 2 files changed, 61 insertions(+), 61 deletions(-) diff --git a/acc/mhd_solver/stencil_defines.h b/acc/mhd_solver/stencil_defines.h index 5124820..ea4c102 100644 --- a/acc/mhd_solver/stencil_defines.h +++ b/acc/mhd_solver/stencil_defines.h @@ -30,8 +30,8 @@ #define LMAGNETIC (1) #define LENTROPY (1) #define LTEMPERATURE (0) -#define LFORCING (1) -#define LUPWD (0) +#define LFORCING (0) +#define LUPWD (1) #define AC_THERMAL_CONDUCTIVITY (AcReal(0.001)) // TODO: make an actual config parameter diff --git a/src/standalone/model/model_rk3.cc b/src/standalone/model/model_rk3.cc index 356f6ad..3471cd1 100644 --- a/src/standalone/model/model_rk3.cc +++ b/src/standalone/model/model_rk3.cc @@ -273,6 +273,53 @@ derzz(const int i, const int j, const int k, const ModelScalar* arr) return second_derivative(pencil, get(AC_inv_dsz)); } +#if LUPWD +static inline ModelScalar +der6x_upwd(const int i, const int j, const int k, const ModelScalar* arr) +{ + ModelScalar inv_ds = get(AC_inv_dsx); + + return ModelScalar(1.0/60.0)*inv_ds* ( + -ModelScalar( 20.0)* arr[IDX(i, j, k)] + +ModelScalar( 15.0)*(arr[IDX(i+1, j, k)] + + arr[IDX(i-1, j, k)]) + -ModelScalar( 6.0)*(arr[IDX(i+2, j, k)] + + arr[IDX(i-2, j, k)]) + + arr[IDX(i+3, j, k)] + + arr[IDX(i-3, j, k)]); +} + +static inline ModelScalar +der6y_upwd(const int i, const int j, const int k, const ModelScalar* arr) +{ + ModelScalar inv_ds = get(AC_inv_dsy); + + return ModelScalar(1.0/60.0)*inv_ds* ( + -ModelScalar( 20.0)* arr[IDX(i, j, k)] + +ModelScalar( 15.0)*(arr[IDX(i, j+1, k)] + + arr[IDX(i, j-1, k)]) + -ModelScalar( 6.0)*(arr[IDX(i, j+2, k)] + + arr[IDX(i, j-2, k)]) + + arr[IDX(i, j+3, k)] + + arr[IDX(i, j-3, k)]); +} + +static inline ModelScalar +der6z_upwd(const int i, const int j, const int k, const ModelScalar* arr) +{ + ModelScalar inv_ds = get(AC_inv_dsz); + + return ModelScalar(1.0/60.0)*inv_ds* ( + -ModelScalar( 20.0)* arr[IDX(i, j, k )] + +ModelScalar( 15.0)*(arr[IDX(i, j, k+1)] + + arr[IDX(i, j, k-1)]) + -ModelScalar( 6.0)*(arr[IDX(i, j, k+2)] + + arr[IDX(i, j, k-2)]) + + arr[IDX(i, j, k+3)] + + arr[IDX(i, j, k-3)]); +} +#endif + static inline ModelScalar compute_value(const int i, const int j, const int k, const ModelScalar* arr) { @@ -354,54 +401,6 @@ gradients(const ModelVectorData& data) return (ModelMatrix){gradient(data.x), gradient(data.y), gradient(data.z)}; } -#if LUPWD - -Preprocessed Scalar -der6x_upwd(in Scalar vertex) -{ - Scalar inv_ds = DCONST_REAL(AC_inv_dsx); - - return (Scalar){ Scalar(1.0/60.0)*inv_ds* ( - - Scalar(20.0)* vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z] - + Scalar(15.0)*(vertex[vertexIdx.x+1, vertexIdx.y, vertexIdx.z] - + vertex[vertexIdx.x-1, vertexIdx.y, vertexIdx.z]) - - Scalar( 6.0)*(vertex[vertexIdx.x+2, vertexIdx.y, vertexIdx.z] - + vertex[vertexIdx.x-2, vertexIdx.y, vertexIdx.z]) - + vertex[vertexIdx.x+3, vertexIdx.y, vertexIdx.z] - + vertex[vertexIdx.x-3, vertexIdx.y, vertexIdx.z])}; -} - -Preprocessed Scalar -der6y_upwd(in Scalar vertex) -{ - Scalar inv_ds = DCONST_REAL(AC_inv_dsy); - - return (Scalar){ Scalar(1.0/60.0)*inv_ds* ( - -Scalar( 20.0)* vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z] - +Scalar( 15.0)*(vertex[vertexIdx.x, vertexIdx.y+1, vertexIdx.z] - + vertex[vertexIdx.x, vertexIdx.y-1, vertexIdx.z]) - -Scalar( 6.0)*(vertex[vertexIdx.x, vertexIdx.y+2, vertexIdx.z] - + vertex[vertexIdx.x, vertexIdx.y-2, vertexIdx.z]) - + vertex[vertexIdx.x, vertexIdx.y+3, vertexIdx.z] - + vertex[vertexIdx.x, vertexIdx.y-3, vertexIdx.z])}; -} - -Preprocessed Scalar -der6z_upwd(in Scalar vertex) -{ - Scalar inv_ds = DCONST_REAL(AC_inv_dsz); - - return (Scalar){ Scalar(1.0/60.0)*inv_ds* ( - -Scalar( 20.0)* vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z] - +Scalar( 15.0)*(vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z+1] - + vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z-1]) - -Scalar( 6.0)*(vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z+2] - + vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z-2]) - + vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z+3] - + vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z-3])}; -} - -#endif /* @@ -500,17 +499,6 @@ curl(const ModelVectorData& vec) gradient(vec.y).x - gradient(vec.x).y}; } -#if LUPWD -Scalar -upwd_der6(in Vector uu, in Scalar lnrho) -{ - Scalar uux = fabs(value(uu).x); - Scalar uuy = fabs(value(uu).y); - Scalar uuz = fabs(value(uu).z); - return (Scalar){uux*der6x_upwd(lnrho) + uuy*der6y_upwd(lnrho) + uuz*der6z_upwd(lnrho)}; -} -#endif - static inline ModelVector gradient_of_divergence(const ModelVectorData& vec) { @@ -563,6 +551,18 @@ contract(const ModelMatrix& mat) * Stencil Processing Stage (equations) * ============================================================================= */ + +#if LUPWD +ModelScalar +upwd_der6(const ModelVectorData& uu, const ModelScalarData& lnrho) +{ + ModelScalar uux = fabs(value(uu).x); + ModelScalar uuy = fabs(value(uu).y); + ModelScalar uuz = fabs(value(uu).z); + return uux*der6x_upwd(lnrho) + uuy*der6y_upwd(lnrho) + uuz*der6z_upwd(lnrho); +} +#endif + static inline ModelScalar continuity(const ModelVectorData& uu, const ModelScalarData& lnrho) {