diff --git a/acc/mhd_solver/stencil_defines.h b/acc/mhd_solver/stencil_defines.h
index b86d287..0ab07c2 100644
--- a/acc/mhd_solver/stencil_defines.h
+++ b/acc/mhd_solver/stencil_defines.h
@@ -1,3 +1,23 @@
+/*
+ Copyright (C) 2014-2019, Johannes Pekkilae, Miikka Vaeisalae.
+
+ This file is part of Astaroth.
+
+ Astaroth is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ Astaroth is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with Astaroth. If not, see .
+*/
+#pragma once
+
/*
* =============================================================================
* Logical switches
@@ -26,8 +46,11 @@
FUNC(AC_max_steps), \
FUNC(AC_save_steps), \
FUNC(AC_bin_steps), \
- FUNC(AC_bc_type)
-#define AC_FOR_REAL_PARAM_TYPES(FUNC)\
+ FUNC(AC_bc_type),
+
+#define AC_FOR_USER_INT3_PARAM_TYPES(FUNC)
+
+#define AC_FOR_USER_REAL_PARAM_TYPES(FUNC)\
/* cparams */\
FUNC(AC_dsx), \
FUNC(AC_dsy), \
@@ -97,7 +120,9 @@
FUNC(AC_cs2_sound), \
FUNC(AC_inv_dsx), \
FUNC(AC_inv_dsy), \
- FUNC(AC_inv_dsz)
+ FUNC(AC_inv_dsz),
+
+#define AC_FOR_USER_REAL3_PARAM_TYPES(FUNC)
// clang-format on
/*
@@ -106,49 +131,33 @@
* =============================================================================
*/
// clang-format off
-#if LDENSITY
-#define AC_FOR_DENSITY_VTXBUF_HANDLES(FUNC) \
- FUNC(VTXBUF_LNRHO),
-#else
-#define AC_FOR_DENSITY_VTXBUF_HANDLES(FUNC)
-#endif
-
-#if LHYDRO
-#define AC_FOR_HYDRO_VTXBUF_HANDLES(FUNC) \
- FUNC(VTXBUF_UUX), \
- FUNC(VTXBUF_UUY), \
- FUNC(VTXBUF_UUZ),
-#else
-#define AC_FOR_HYDRO_VTXBUF_HANDLES(FUNC)
-#endif
-
-#if LMAGNETIC
-#define AC_FOR_MAGNETIC_VTXBUF_HANDLES(FUNC) \
- FUNC(VTXBUF_AX), \
- FUNC(VTXBUF_AY), \
- FUNC(VTXBUF_AZ),
-#else
-#define AC_FOR_MAGNETIC_VTXBUF_HANDLES(FUNC)
-#endif
-
#if LENTROPY
-#define AC_FOR_ENTROPY_VTXBUF_HANDLES(FUNC) \
- FUNC(VTXBUF_ENTROPY),
+#define AC_FOR_VTXBUF_HANDLES(FUNC) \
+ FUNC(VTXBUF_LNRHO), \
+ FUNC(VTXBUF_UUX), \
+ FUNC(VTXBUF_UUY), \
+ FUNC(VTXBUF_UUZ), \
+ FUNC(VTXBUF_AX), \
+ FUNC(VTXBUF_AY), \
+ FUNC(VTXBUF_AZ), \
+ FUNC(VTXBUF_ENTROPY),
+#elif LMAGNETIC
+#define AC_FOR_VTXBUF_HANDLES(FUNC) \
+ FUNC(VTXBUF_LNRHO), \
+ FUNC(VTXBUF_UUX), \
+ FUNC(VTXBUF_UUY), \
+ FUNC(VTXBUF_UUZ), \
+ FUNC(VTXBUF_AX), \
+ FUNC(VTXBUF_AY), \
+ FUNC(VTXBUF_AZ),
+#elif LHYDRO
+#define AC_FOR_VTXBUF_HANDLES(FUNC) \
+ FUNC(VTXBUF_LNRHO), \
+ FUNC(VTXBUF_UUX), \
+ FUNC(VTXBUF_UUY), \
+ FUNC(VTXBUF_UUZ),
#else
-#define AC_FOR_ENTROPY_VTXBUF_HANDLES(FUNC)
+#define AC_FOR_VTXBUF_HANDLES(FUNC) \
+ FUNC(VTXBUF_LNRHO),
#endif
-
-//MR: Temperature must not have an additional variable slot, but should sit on the
-// same as entropy.
-#if LTEMPERATURE
- #define AC_FOR_TEMPERATURE_VTXBUF_HANDLES(FUNC)\
- FUNC(VTXBUF_TEMPERATURE),
-#else
- #define AC_FOR_TEMPERATURE_VTXBUF_HANDLES(FUNC)
-#endif
-
-#define AC_FOR_VTXBUF_HANDLES(FUNC) AC_FOR_HYDRO_VTXBUF_HANDLES(FUNC) \
- AC_FOR_DENSITY_VTXBUF_HANDLES(FUNC) \
- AC_FOR_ENTROPY_VTXBUF_HANDLES(FUNC) \
- AC_FOR_MAGNETIC_VTXBUF_HANDLES(FUNC) \
// clang-format on
diff --git a/acc/samples/common_header.h b/acc/samples/common_header.h
index 5e4e4e6..4701348 100644
--- a/acc/samples/common_header.h
+++ b/acc/samples/common_header.h
@@ -273,20 +273,20 @@ typedef enum {
*/
typedef enum {
AC_FOR_INT_PARAM_TYPES(AC_GEN_ID),
- NUM_INT_PARAM_TYPES
+ NUM_INT_PARAMS
} AcIntParam;
typedef enum {
AC_FOR_REAL_PARAM_TYPES(AC_GEN_ID),
- NUM_REAL_PARAM_TYPES
+ NUM_REAL_PARAMS
} AcRealParam;
extern const char* intparam_names[]; // Defined in astaroth.cu
extern const char* realparam_names[]; // Defined in astaroth.cu
typedef struct {
- int int_params[NUM_INT_PARAM_TYPES];
- AcReal real_params[NUM_REAL_PARAM_TYPES];
+ int int_params[NUM_INT_PARAMS];
+ AcReal real_params[NUM_REAL_PARAMS];
} AcMeshInfo;
/*
@@ -316,21 +316,21 @@ typedef struct {
AcMeshInfo info;
} AcMesh;
-#define AC_VTXBUF_SIZE(mesh_info) \
+#define acVertexBufferSize(mesh_info) \
((size_t)(mesh_info.int_params[AC_mx] * mesh_info.int_params[AC_my] * \
mesh_info.int_params[AC_mz]))
-#define AC_VTXBUF_SIZE_BYTES(mesh_info) \
- (sizeof(AcReal) * AC_VTXBUF_SIZE(mesh_info))
+#define acVertexBufferSizeBytes(mesh_info) \
+ (sizeof(AcReal) * acVertexBufferSize(mesh_info))
-#define AC_VTXBUF_COMPDOMAIN_SIZE(mesh_info) \
+#define acVertexBufferCompdomainSize(mesh_info) \
(mesh_info.int_params[AC_nx] * mesh_info.int_params[AC_ny] * \
mesh_info.int_params[AC_nz])
-#define AC_VTXBUF_COMPDOMAIN_SIZE_BYTES(mesh_info) \
- (sizeof(AcReal) * AC_VTXBUF_COMPDOMAIN_SIZE(mesh_info))
+#define acVertexBufferCompdomainSizeBytes(mesh_info) \
+ (sizeof(AcReal) * acVertexBufferCompdomainSize(mesh_info))
-#define AC_VTXBUF_IDX(i, j, k, mesh_info) \
+#define acVertexBufferIdx(i, j, k, mesh_info) \
((i) + (j)*mesh_info.int_params[AC_mx] + \
(k)*mesh_info.int_params[AC_mx] * mesh_info.int_params[AC_my])
diff --git a/doc/manual/manual.md b/doc/manual/manual.md
index 007b1bc..97eac88 100644
--- a/doc/manual/manual.md
+++ b/doc/manual/manual.md
@@ -56,7 +56,7 @@ Saving output binaries is not enabled yet.
- All essential tructs, macros and enumerators are found in astaroth.h for better reference.
-- In the case there is changes in the data layout, better use macro `AC_VTXBUF_IDX(i, j, k, mesh_info)`which transform indices from 3D to 1D. Therefore no need to start writing `i + j * mesh_info.int_params[AC_mx] + ...` which would affect the code readability.
+- In the case there is changes in the data layout, better use macro `acVertexBufferIdx(i, j, k, mesh_info)`which transform indices from 3D to 1D. Therefore no need to start writing `i + j * mesh_info.int_params[AC_mx] + ...` which would affect the code readability.
- AcReal on generic floating point real number type used everywhere in the code. Currently can be either `float` or `double`. Possibly in the future also `half` or `long double` could become available.
@@ -96,7 +96,7 @@ printf("nx: %d, dsx %f\n",
mesh->info.int_params[AC_nx],
double(mesh->info.real_params[AC_dsx]));
printf("First vertex of the computational domain: %f\n",
-double(mesh->vertex_buffer[VTXBUF_LNRHO][AC_VTXBUF_IDX(3, 3, 3, mesh_info)]));
+double(mesh->vertex_buffer[VTXBUF_LNRHO][acVertexBufferIdx(3, 3, 3, mesh_info)]));
```
diff --git a/include/astaroth.h b/include/astaroth.h
index 5c53238..be0081e 100644
--- a/include/astaroth.h
+++ b/include/astaroth.h
@@ -16,209 +16,13 @@
You should have received a copy of the GNU General Public License
along with Astaroth. If not, see .
*/
-
-/**
- * @file
- * \brief Brief info.
- *
- * Provides an interface to Astaroth. Contains all the necessary configuration
- * structs and functions for running the code on multiple GPUs.
- *
- * All interface functions declared here (such as acInit()) operate all GPUs
- * available in the node under the hood, and the user does not need any
- * information about the decomposition, synchronization or such to use these
- * functions.
- *
- */
#pragma once
-/* Prevent name mangling */
#ifdef __cplusplus
extern "C" {
#endif
-#include // FLT_EPSILON, etc
-#include // size_t
-#include // CUDA vector types (float4, etc)
-
-/*
- * =============================================================================
- * Flags for auto-optimization
- * =============================================================================
- */
-#define AUTO_OPTIMIZE (0) // DEPRECATED TODO remove
-#define BOUNDCONDS_OPTIMIZE (0)
-#define GENERATE_BENCHMARK_DATA (0)
-#define VERBOSE_PRINTING (1)
-
-// Device info
-#define REGISTERS_PER_THREAD (255)
-#define MAX_REGISTERS_PER_BLOCK (65536)
-#define MAX_THREADS_PER_BLOCK (1024)
-#define WARP_SIZE (32)
-/*
- * =============================================================================
- * Compile-time constants used during simulation (user definable)
- * =============================================================================
- */
-// USER_PROVIDED_DEFINES must be defined in user.h if the user wants to override the following
-// logical switches
-#include "user.h"
-
-// clang-format off
-#ifndef USER_PROVIDED_DEFINES
- #include "stencil_defines.h"
-#endif
-// clang-format on
-
-/*
- * =============================================================================
- * Built-in parameters
- * =============================================================================
- */
-// clang-format off
-#define AC_FOR_BUILTIN_INT_PARAM_TYPES(FUNC)\
- /* cparams */\
- FUNC(AC_nx), \
- FUNC(AC_ny), \
- FUNC(AC_nz), \
- FUNC(AC_mx), \
- FUNC(AC_my), \
- FUNC(AC_mz), \
- FUNC(AC_nx_min), \
- FUNC(AC_ny_min), \
- FUNC(AC_nz_min), \
- FUNC(AC_nx_max), \
- FUNC(AC_ny_max), \
- FUNC(AC_nz_max), \
- /* Additional */\
- FUNC(AC_mxy),\
- FUNC(AC_nxy),\
- FUNC(AC_nxyz),
-// clang-format on
-
-/*
- * =============================================================================
- * Single/double precision switch
- * =============================================================================
- */
-// clang-format off
-#if AC_DOUBLE_PRECISION == 1
- typedef double AcReal;
- typedef double3 AcReal3;
- #define AC_REAL_MAX (DBL_MAX)
- #define AC_REAL_MIN (DBL_MIN)
- #define AC_REAL_EPSILON (DBL_EPSILON)
-#else
- typedef float AcReal;
- typedef float3 AcReal3;
- #define AC_REAL_MAX (FLT_MAX)
- #define AC_REAL_MIN (FLT_MIN)
- #define AC_REAL_EPSILON (FLT_EPSILON)
-#endif
-// clang-format on
-
-typedef struct {
- AcReal3 row[3];
-} AcMatrix;
-
-/*
- * =============================================================================
- * Helper macros
- * =============================================================================
- */
-#define AC_GEN_ID(X) X
-#define AC_GEN_STR(X) #X
-
-/*
- * =============================================================================
- * Error codes
- * =============================================================================
- */
-typedef enum { AC_SUCCESS = 0, AC_FAILURE = 1 } AcResult;
-
-/*
- * =============================================================================
- * Reduction types
- * =============================================================================
- */
-typedef enum { RTYPE_MAX, RTYPE_MIN, RTYPE_RMS, RTYPE_RMS_EXP, NUM_REDUCTION_TYPES } ReductionType;
-
-/*
- * =============================================================================
- * Definitions for the enums and structs for AcMeshInfo (DO NOT TOUCH)
- * =============================================================================
- */
-
-typedef enum {
- AC_FOR_BUILTIN_INT_PARAM_TYPES(AC_GEN_ID) //
- AC_FOR_USER_INT_PARAM_TYPES(AC_GEN_ID), //
- NUM_INT_PARAM_TYPES
-} AcIntParam;
-
-typedef enum { AC_FOR_REAL_PARAM_TYPES(AC_GEN_ID), NUM_REAL_PARAM_TYPES } AcRealParam;
-// typedef enum { AC_FOR_VEC_PARAM_TYPES(AC_GEN_ID), NUM_VEC_PARAM_TYPES } AcVecParam;
-
-extern const char* intparam_names[]; // Defined in astaroth.cu
-extern const char* realparam_names[]; // Defined in astaroth.cu
-
-typedef struct {
- int int_params[NUM_INT_PARAM_TYPES];
- AcReal real_params[NUM_REAL_PARAM_TYPES];
- // AcReal* vec_params[NUM_VEC_PARAM_TYPES];
-} AcMeshInfo;
-
-/*
- * =============================================================================
- * Definitions for the enums and structs for AcMesh (DO NOT TOUCH)
- * =============================================================================
- */
-typedef enum { AC_FOR_VTXBUF_HANDLES(AC_GEN_ID) NUM_VTXBUF_HANDLES } VertexBufferHandle;
-
-extern const char* vtxbuf_names[]; // Defined in astaroth.cu
-
-/*
-typedef struct {
- AcReal* data;
-} VertexBuffer;
-*/
-
-// NOTE: there's no particular benefit declaring AcMesh a class, since
-// a library user may already have allocated memory for the vertex_buffers.
-// But then we would allocate memory again when the user wants to start
-// filling the class with data. => Its better to consider AcMesh as a
-// payload-only struct
-typedef struct {
- AcReal* vertex_buffer[NUM_VTXBUF_HANDLES];
- AcMeshInfo info;
-} AcMesh;
-
-#define AC_VTXBUF_SIZE(mesh_info) \
- ((size_t)(mesh_info.int_params[AC_mx] * mesh_info.int_params[AC_my] * \
- mesh_info.int_params[AC_mz]))
-
-#define AC_VTXBUF_SIZE_BYTES(mesh_info) (sizeof(AcReal) * AC_VTXBUF_SIZE(mesh_info))
-
-#define AC_VTXBUF_COMPDOMAIN_SIZE(mesh_info) \
- (mesh_info.int_params[AC_nx] * mesh_info.int_params[AC_ny] * mesh_info.int_params[AC_nz])
-
-#define AC_VTXBUF_COMPDOMAIN_SIZE_BYTES(mesh_info) \
- (sizeof(AcReal) * AC_VTXBUF_COMPDOMAIN_SIZE(mesh_info))
-
-#define AC_VTXBUF_IDX(i, j, k, mesh_info) \
- ((i) + (j)*mesh_info.int_params[AC_mx] + \
- (k)*mesh_info.int_params[AC_mx] * mesh_info.int_params[AC_my])
-
-/*
- * =============================================================================
- * Astaroth interface: Basic functions. Synchronous.
- * =============================================================================
- */
-typedef enum {
- STREAM_DEFAULT,
- NUM_STREAM_TYPES, //
- STREAM_ALL
-} StreamType;
+#include "astaroth_defines.h"
/** Checks whether there are any CUDA devices available. Returns AC_SUCCESS if there is 1 or more,
* AC_FAILURE otherwise. */
@@ -304,7 +108,6 @@ AcResult acIntegrateStepWithOffsetAsync(const int& isubstep, const AcReal& dt, c
AcResult acBoundcondStep(void);
AcResult acBoundcondStepAsync(const StreamType stream);
-/* End extern "C" */
#ifdef __cplusplus
-}
+} // extern "C"
#endif
diff --git a/include/astaroth_defines.h b/include/astaroth_defines.h
new file mode 100644
index 0000000..ac9804b
--- /dev/null
+++ b/include/astaroth_defines.h
@@ -0,0 +1,226 @@
+/*
+ Copyright (C) 2014-2019, Johannes Pekkilae, Miikka Vaeisalae.
+
+ This file is part of Astaroth.
+
+ Astaroth is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ Astaroth is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with Astaroth. If not, see .
+*/
+#pragma once
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+#include // FLT_EPSILON, etc
+#include // size_t
+#include // CUDA vector types (float4, etc)
+
+#include "stencil_defines.h"
+
+// Library flags
+#define VERBOSE_PRINTING (1)
+
+// Built-in types and parameters
+#if AC_DOUBLE_PRECISION == 1
+typedef double AcReal;
+typedef double3 AcReal3;
+#define AC_REAL_MAX (DBL_MAX)
+#define AC_REAL_MIN (DBL_MIN)
+#define AC_REAL_EPSILON (DBL_EPSILON)
+#else
+typedef float AcReal;
+typedef float3 AcReal3;
+#define AC_REAL_MAX (FLT_MAX)
+#define AC_REAL_MIN (FLT_MIN)
+#define AC_REAL_EPSILON (FLT_EPSILON)
+#endif
+
+typedef struct {
+ AcReal3 row[3];
+} AcMatrix;
+
+// clang-format off
+#define AC_FOR_BUILTIN_INT_PARAM_TYPES(FUNC)\
+ FUNC(AC_nx), \
+ FUNC(AC_ny), \
+ FUNC(AC_nz), \
+ FUNC(AC_mx), \
+ FUNC(AC_my), \
+ FUNC(AC_mz), \
+ FUNC(AC_nx_min), \
+ FUNC(AC_ny_min), \
+ FUNC(AC_nz_min), \
+ FUNC(AC_nx_max), \
+ FUNC(AC_ny_max), \
+ FUNC(AC_nz_max), \
+ FUNC(AC_mxy),\
+ FUNC(AC_nxy),\
+ FUNC(AC_nxyz),\
+
+#define AC_FOR_BUILTIN_INT3_PARAM_TYPES(FUNC)
+
+#define AC_FOR_BUILTIN_REAL_PARAM_TYPES(FUNC)
+
+#define AC_FOR_BUILTIN_REAL3_PARAM_TYPES(FUNC)
+// clang-format on
+
+typedef enum { AC_SUCCESS = 0, AC_FAILURE = 1 } AcResult;
+
+typedef enum { RTYPE_MAX, RTYPE_MIN, RTYPE_RMS, RTYPE_RMS_EXP, NUM_REDUCTION_TYPES } ReductionType;
+
+typedef enum {
+ STREAM_DEFAULT,
+ NUM_STREAM_TYPES, //
+ STREAM_ALL
+} StreamType;
+
+#define AC_GEN_ID(X) X
+typedef enum {
+ AC_FOR_BUILTIN_INT_PARAM_TYPES(AC_GEN_ID) //
+ AC_FOR_USER_INT_PARAM_TYPES(AC_GEN_ID) //
+ NUM_INT_PARAMS
+} AcIntParam;
+
+typedef enum {
+ AC_FOR_BUILTIN_INT3_PARAM_TYPES(AC_GEN_ID) //
+ AC_FOR_USER_INT3_PARAM_TYPES(AC_GEN_ID) //
+ NUM_INT3_PARAMS
+} AcInt3Param;
+
+typedef enum {
+ AC_FOR_BUILTIN_REAL_PARAM_TYPES(AC_GEN_ID) //
+ AC_FOR_USER_REAL_PARAM_TYPES(AC_GEN_ID) //
+ NUM_REAL_PARAMS
+} AcRealParam;
+
+typedef enum {
+ AC_FOR_BUILTIN_REAL3_PARAM_TYPES(AC_GEN_ID) //
+ AC_FOR_USER_REAL3_PARAM_TYPES(AC_GEN_ID) //
+ NUM_REAL3_PARAMS
+} AcReal3Param;
+
+typedef enum {
+ AC_FOR_VTXBUF_HANDLES(AC_GEN_ID) //
+ NUM_VTXBUF_HANDLES
+} VertexBufferHandle;
+#undef AC_GEN_ID
+
+extern const char* intparam_names[];
+extern const char* int3param_names[];
+extern const char* realparam_names[];
+extern const char* real3param_names[];
+extern const char* vtxbuf_names[];
+
+typedef struct {
+ int int_params[NUM_INT_PARAMS];
+ int3 int3_params[NUM_INT3_PARAMS];
+ AcReal real_params[NUM_REAL_PARAMS];
+ AcReal3 real3_params[NUM_REAL3_PARAMS];
+} AcMeshInfo;
+
+typedef struct {
+ AcReal* vertex_buffer[NUM_VTXBUF_HANDLES];
+ AcMeshInfo info;
+} AcMesh;
+
+/*
+ * =============================================================================
+ * Helper functions
+ * =============================================================================
+ */
+static inline size_t
+acVertexBufferSize(const AcMeshInfo& info)
+{
+ return info.int_params[AC_mx] * info.int_params[AC_my] * info.int_params[AC_mz];
+}
+
+static inline size_t
+acVertexBufferSizeBytes(const AcMeshInfo& info)
+{
+ return sizeof(AcReal) * acVertexBufferSize(info);
+}
+
+static inline size_t
+acVertexBufferCompdomainSize(const AcMeshInfo& info)
+{
+ return info.int_params[AC_nx] * info.int_params[AC_ny] * info.int_params[AC_nz];
+}
+
+static inline size_t
+acVertexBufferCompdomainSizeBytes(const AcMeshInfo& info)
+{
+ return sizeof(AcReal) * acVertexBufferCompdomainSize(info);
+}
+
+static inline size_t
+acVertexBufferIdx(const int i, const int j, const int k, const AcMeshInfo& info)
+{
+ return i + //
+ j * info.int_params[AC_mx] + //
+ k * info.int_params[AC_mx] * info.int_params[AC_my];
+}
+
+/*
+static inline int
+acGetParam(const AcMeshInfo& info, const AcIntParam param)
+{
+ return info.int_params[param];
+}
+
+static inline int3
+acGetParam(const AcMeshInfo& info, const AcInt3Param param)
+{
+ return info.int3_params[param];
+}
+
+static inline AcReal
+acGetParam(const AcMeshInfo& info, const AcRealParam param)
+{
+ return info.real_params[param];
+}
+
+static inline AcReal3
+acGetParam(const AcMeshInfo& info, const AcReal3Param param)
+{
+ return info.real3_params[param];
+}
+
+static inline void
+acSetParam(const AcIntParam param, const int value, AcMeshInfo* info)
+{
+ info->int_params[param] = value;
+}
+
+static inline void
+acSetParam(const AcInt3Param param, const int3 value, AcMeshInfo* info)
+{
+ info->int3_params[param] = value;
+}
+
+static inline void
+acSetParam(const AcRealParam param, const AcReal value, AcMeshInfo* info)
+{
+ info->real_params[param] = value;
+}
+
+static inline void
+acSetParam(const AcReal3Param param, const AcReal3 value, AcMeshInfo* info)
+{
+ info->real3_params[param] = value;
+}
+*/
+
+#ifdef __cplusplus
+} // extern "C"
+#endif
diff --git a/src/core/astaroth.cu b/src/core/astaroth.cu
index 38b2863..fadf4c0 100644
--- a/src/core/astaroth.cu
+++ b/src/core/astaroth.cu
@@ -130,10 +130,17 @@
#include "math_utils.h" // sum for reductions
#include "standalone/config_loader.h" // update_config
-const char* intparam_names[] = {AC_FOR_BUILTIN_INT_PARAM_TYPES(AC_GEN_STR)
- AC_FOR_USER_INT_PARAM_TYPES(AC_GEN_STR)};
-const char* realparam_names[] = {AC_FOR_REAL_PARAM_TYPES(AC_GEN_STR)};
-const char* vtxbuf_names[] = {AC_FOR_VTXBUF_HANDLES(AC_GEN_STR)};
+#define AC_GEN_STR(X) #X
+const char* intparam_names[] = {AC_FOR_BUILTIN_INT_PARAM_TYPES(AC_GEN_STR) //
+ AC_FOR_USER_INT_PARAM_TYPES(AC_GEN_STR)};
+const char* int3param_names[] = {AC_FOR_BUILTIN_INT3_PARAM_TYPES(AC_GEN_STR) //
+ AC_FOR_USER_INT3_PARAM_TYPES(AC_GEN_STR)};
+const char* realparam_names[] = {AC_FOR_BUILTIN_REAL_PARAM_TYPES(AC_GEN_STR) //
+ AC_FOR_USER_REAL_PARAM_TYPES(AC_GEN_STR)};
+const char* real3param_names[] = {AC_FOR_BUILTIN_REAL3_PARAM_TYPES(AC_GEN_STR) //
+ AC_FOR_USER_REAL3_PARAM_TYPES(AC_GEN_STR)};
+const char* vtxbuf_names[] = {AC_FOR_VTXBUF_HANDLES(AC_GEN_STR)};
+#undef AC_GEN_STR
static const int MAX_NUM_DEVICES = 32;
static int num_devices = 0;
@@ -559,7 +566,7 @@ acLoadWithOffset(const AcMesh& host_mesh, const int3& src, const int num_vertice
AcResult
acLoad(const AcMesh& host_mesh)
{
- acLoadWithOffset(host_mesh, (int3){0, 0, 0}, AC_VTXBUF_SIZE(host_mesh.info));
+ acLoadWithOffset(host_mesh, (int3){0, 0, 0}, acVertexBufferSize(host_mesh.info));
acSynchronizeStream(STREAM_ALL);
return AC_SUCCESS;
}
@@ -598,7 +605,7 @@ acStoreWithOffset(const int3& src, const int num_vertices, AcMesh* host_mesh)
AcResult
acStore(AcMesh* host_mesh)
{
- acStoreWithOffset((int3){0, 0, 0}, AC_VTXBUF_SIZE(host_mesh->info), host_mesh);
+ acStoreWithOffset((int3){0, 0, 0}, acVertexBufferSize(host_mesh->info), host_mesh);
acSynchronizeStream(STREAM_ALL);
return AC_SUCCESS;
}
diff --git a/src/core/device.cu b/src/core/device.cu
index d86311b..7b624e3 100644
--- a/src/core/device.cu
+++ b/src/core/device.cu
@@ -28,6 +28,12 @@
#include "errchk.h"
+// Device info
+#define REGISTERS_PER_THREAD (255)
+#define MAX_REGISTERS_PER_BLOCK (65536)
+#define MAX_THREADS_PER_BLOCK (1024)
+#define WARP_SIZE (32)
+
typedef struct {
AcReal* in[NUM_VTXBUF_HANDLES];
AcReal* out[NUM_VTXBUF_HANDLES];
@@ -159,13 +165,13 @@ createDevice(const int id, const AcMeshInfo device_config, Device* device_handle
}
// Memory
- const size_t vba_size_bytes = AC_VTXBUF_SIZE_BYTES(device_config);
+ const size_t vba_size_bytes = acVertexBufferSizeBytes(device_config);
for (int i = 0; i < NUM_VTXBUF_HANDLES; ++i) {
ERRCHK_CUDA_ALWAYS(cudaMalloc(&device->vba.in[i], vba_size_bytes));
ERRCHK_CUDA_ALWAYS(cudaMalloc(&device->vba.out[i], vba_size_bytes));
}
ERRCHK_CUDA_ALWAYS(
- cudaMalloc(&device->reduce_scratchpad, AC_VTXBUF_COMPDOMAIN_SIZE_BYTES(device_config)));
+ cudaMalloc(&device->reduce_scratchpad, acVertexBufferCompdomainSizeBytes(device_config)));
ERRCHK_CUDA_ALWAYS(cudaMalloc(&device->reduce_result, sizeof(AcReal)));
#if PACKED_DATA_TRANSFERS
@@ -335,8 +341,8 @@ AcResult
copyMeshToDevice(const Device device, const StreamType stream_type, const AcMesh& host_mesh,
const int3& src, const int3& dst, const int num_vertices)
{
- const size_t src_idx = AC_VTXBUF_IDX(src.x, src.y, src.z, host_mesh.info);
- const size_t dst_idx = AC_VTXBUF_IDX(dst.x, dst.y, dst.z, device->local_config);
+ const size_t src_idx = acVertexBufferIdx(src.x, src.y, src.z, host_mesh.info);
+ const size_t dst_idx = acVertexBufferIdx(dst.x, dst.y, dst.z, device->local_config);
for (int i = 0; i < NUM_VTXBUF_HANDLES; ++i) {
loadWithOffset(device, stream_type, &host_mesh.vertex_buffer[i][src_idx],
num_vertices * sizeof(AcReal), &device->vba.in[i][dst_idx]);
@@ -348,8 +354,8 @@ AcResult
copyMeshToHost(const Device device, const StreamType stream_type, const int3& src, const int3& dst,
const int num_vertices, AcMesh* host_mesh)
{
- const size_t src_idx = AC_VTXBUF_IDX(src.x, src.y, src.z, device->local_config);
- const size_t dst_idx = AC_VTXBUF_IDX(dst.x, dst.y, dst.z, host_mesh->info);
+ const size_t src_idx = acVertexBufferIdx(src.x, src.y, src.z, device->local_config);
+ const size_t dst_idx = acVertexBufferIdx(dst.x, dst.y, dst.z, host_mesh->info);
for (int i = 0; i < NUM_VTXBUF_HANDLES; ++i) {
storeWithOffset(device, stream_type, &device->vba.in[i][src_idx],
num_vertices * sizeof(AcReal), &host_mesh->vertex_buffer[i][dst_idx]);
@@ -362,8 +368,8 @@ copyMeshDeviceToDevice(const Device src_device, const StreamType stream_type, co
Device dst_device, const int3& dst, const int num_vertices)
{
cudaSetDevice(src_device->id);
- const size_t src_idx = AC_VTXBUF_IDX(src.x, src.y, src.z, src_device->local_config);
- const size_t dst_idx = AC_VTXBUF_IDX(dst.x, dst.y, dst.z, dst_device->local_config);
+ const size_t src_idx = acVertexBufferIdx(src.x, src.y, src.z, src_device->local_config);
+ const size_t dst_idx = acVertexBufferIdx(dst.x, dst.y, dst.z, dst_device->local_config);
for (int i = 0; i < NUM_VTXBUF_HANDLES; ++i) {
ERRCHK_CUDA(cudaMemcpyPeerAsync(&dst_device->vba.in[i][dst_idx], dst_device->id,
diff --git a/src/standalone/autotest.cc b/src/standalone/autotest.cc
index bb8b55f..c7aec10 100644
--- a/src/standalone/autotest.cc
+++ b/src/standalone/autotest.cc
@@ -318,7 +318,7 @@ verify_meshes(const ModelMesh& model, const AcMesh& candidate)
const ModelScalar range = get_data_range(model);
for (int w = 0; w < NUM_VTXBUF_HANDLES; ++w) {
- const size_t n = AC_VTXBUF_SIZE(model.info);
+ const size_t n = acVertexBufferSize(model.info);
// Maximum errors
ErrorInfo max_abs_error = ErrorInfo();
@@ -555,7 +555,7 @@ get_max_abs_error_mesh(const ModelMesh& model_mesh, const AcMesh& candidate_mesh
error.abs_error = -1;
for (size_t j = 0; j < NUM_VTXBUF_HANDLES; ++j) {
- for (size_t i = 0; i < AC_VTXBUF_SIZE(model_mesh.info); ++i) {
+ for (size_t i = 0; i < acVertexBufferSize(model_mesh.info); ++i) {
Error curr_error = get_error(model_mesh.vertex_buffer[j][i],
candidate_mesh.vertex_buffer[j][i]);
if (curr_error.abs_error > error.abs_error)
@@ -574,7 +574,7 @@ get_maximum_magnitude(const ModelScalar* field, const AcMeshInfo info)
{
ModelScalar maximum = -INFINITY;
- for (size_t i = 0; i < AC_VTXBUF_SIZE(info); ++i)
+ for (size_t i = 0; i < acVertexBufferSize(info); ++i)
maximum = max(maximum, fabsl(field[i]));
return maximum;
@@ -585,7 +585,7 @@ get_minimum_magnitude(const ModelScalar* field, const AcMeshInfo info)
{
ModelScalar minimum = INFINITY;
- for (size_t i = 0; i < AC_VTXBUF_SIZE(info); ++i)
+ for (size_t i = 0; i < acVertexBufferSize(info); ++i)
minimum = min(minimum, fabsl(field[i]));
return minimum;
@@ -601,7 +601,7 @@ get_max_abs_error_vtxbuf(const VertexBufferHandle vtxbuf_handle, const ModelMesh
Error error;
error.abs_error = -1;
- for (size_t i = 0; i < AC_VTXBUF_SIZE(model_mesh.info); ++i) {
+ for (size_t i = 0; i < acVertexBufferSize(model_mesh.info); ++i) {
Error curr_error = get_error(model_vtxbuf[i], candidate_vtxbuf[i]);
diff --git a/src/standalone/config_loader.cc b/src/standalone/config_loader.cc
index 3eb33a8..7a0a509 100644
--- a/src/standalone/config_loader.cc
+++ b/src/standalone/config_loader.cc
@@ -37,9 +37,9 @@
static inline void
print(const AcMeshInfo& config)
{
- for (int i = 0; i < NUM_INT_PARAM_TYPES; ++i)
+ for (int i = 0; i < NUM_INT_PARAMS; ++i)
printf("[%s]: %d\n", intparam_names[i], config.int_params[i]);
- for (int i = 0; i < NUM_REAL_PARAM_TYPES; ++i)
+ for (int i = 0; i < NUM_REAL_PARAMS; ++i)
printf("[%s]: %g\n", realparam_names[i], double(config.real_params[i]));
}
@@ -77,9 +77,9 @@ parse_config(const char* path, AcMeshInfo* config)
continue;
int idx = -1;
- if ((idx = find_str(keyword, intparam_names, NUM_INT_PARAM_TYPES)) >= 0)
+ if ((idx = find_str(keyword, intparam_names, NUM_INT_PARAMS)) >= 0)
config->int_params[idx] = atoi(value);
- else if ((idx = find_str(keyword, realparam_names, NUM_REAL_PARAM_TYPES)) >= 0)
+ else if ((idx = find_str(keyword, realparam_names, NUM_REAL_PARAMS)) >= 0)
config->real_params[idx] = AcReal(atof(value));
}
diff --git a/src/standalone/model/host_memory.cc b/src/standalone/model/host_memory.cc
index 0762b30..5cda923 100644
--- a/src/standalone/model/host_memory.cc
+++ b/src/standalone/model/host_memory.cc
@@ -30,7 +30,9 @@
#include "core/errchk.h"
+#define AC_GEN_STR(X) #X
const char* init_type_names[] = {AC_FOR_INIT_TYPES(AC_GEN_STR)};
+#undef AC_GEN_STR
#define XORIG (AcReal(.5) * mesh->info.int_params[AC_nx] * mesh->info.real_params[AC_dsx])
#define YORIG (AcReal(.5) * mesh->info.int_params[AC_ny] * mesh->info.real_params[AC_dsy])
@@ -43,14 +45,14 @@ static uint64_t ac_rand_next = 1;
static int32_t
ac_rand(void)
{
- ac_rand_next = ac_rand_next * 1103515245 + 12345;
- return (uint32_t)(ac_rand_next/65536) % 32768;
+ ac_rand_next = ac_rand_next * 1103515245 + 12345;
+ return (uint32_t)(ac_rand_next/65536) % 32768;
}
static void
ac_srand(const uint32_t seed)
{
- ac_rand_next = seed;
+ ac_rand_next = seed;
}
*/
@@ -60,7 +62,7 @@ acmesh_create(const AcMeshInfo& mesh_info)
AcMesh* mesh = (AcMesh*)malloc(sizeof(*mesh));
mesh->info = mesh_info;
- const size_t bytes = AC_VTXBUF_SIZE_BYTES(mesh->info);
+ const size_t bytes = acVertexBufferSizeBytes(mesh->info);
for (int i = 0; i < NUM_VTXBUF_HANDLES; ++i) {
mesh->vertex_buffer[VertexBufferHandle(i)] = (AcReal*)malloc(bytes);
ERRCHK(mesh->vertex_buffer[VertexBufferHandle(i)] != NULL);
@@ -70,15 +72,13 @@ acmesh_create(const AcMeshInfo& mesh_info)
}
static void
-vertex_buffer_set(const VertexBufferHandle& key, const AcReal& val,
- AcMesh* mesh)
+vertex_buffer_set(const VertexBufferHandle& key, const AcReal& val, AcMesh* mesh)
{
- const int n = AC_VTXBUF_SIZE(mesh->info);
+ const int n = acVertexBufferSize(mesh->info);
for (int i = 0; i < n; ++i)
mesh->vertex_buffer[key][i] = val;
}
-
/** Inits all fields to 1. Setting the mesh to zero is problematic because some fields are supposed
to be > 0 and the results would vary widely, which leads to loss of precision in the
computations */
@@ -95,13 +95,12 @@ randr(void)
return AcReal(rand()) / AcReal(RAND_MAX);
}
-
void
lnrho_step(AcMesh* mesh)
{
- const int mx = mesh->info.int_params[AC_mx];
- const int my = mesh->info.int_params[AC_my];
- const int mz = mesh->info.int_params[AC_mz];
+ const int mx = mesh->info.int_params[AC_mx];
+ const int my = mesh->info.int_params[AC_my];
+ const int mz = mesh->info.int_params[AC_mz];
// const int nx_min = mesh->info.int_params[AC_nx_min];
// const int nx_max = mesh->info.int_params[AC_nx_max];
@@ -116,22 +115,23 @@ lnrho_step(AcMesh* mesh)
// const AcReal xmax = DX * (nx_max - nx_min) ;
// const AcReal zmax = DZ * (nz_max - nz_min) ;
- // const AcReal lnrho1 = (AcReal) -1.0; // TODO mesh->info.real_params[AC_lnrho1];
- const AcReal lnrho2 = (AcReal) 0.0; // TODO mesh->info.real_params[AC_lnrho2];
- // const AcReal rho1 = (AcReal) exp(lnrho1);
+ // const AcReal lnrho1 = (AcReal) -1.0; // TODO mesh->info.real_params[AC_lnrho1];
+ const AcReal lnrho2 = (AcReal)0.0; // TODO mesh->info.real_params[AC_lnrho2];
+ // const AcReal rho1 = (AcReal) exp(lnrho1);
// const AcReal rho2 = (AcReal) exp(lnrho2);
- // const AcReal k_pert = (AcReal) 1.0; //mesh->info.real_params[AC_k_pert]; //Wamenumber of the perturbation
- // const AcReal k_pert = 4.0; //mesh->info.real_params[AC_k_pert]; //Wamenumber of the perturbation
- //const AcReal ampl_pert = xmax/10.0; // xmax/mesh->info.real_params[AC_pert]; //Amplitude of the perturbation
- // const AcReal ampl_pert = (AcReal) 0.0;//xmax/20.0; // xmax/mesh->info.real_params[AC_pert]; //Amplitude of the perturbation
- // const AcReal two_pi = (AcReal) 6.28318531;
+ // const AcReal k_pert = (AcReal) 1.0; //mesh->info.real_params[AC_k_pert]; //Wamenumber of
+ // the perturbation const AcReal k_pert = 4.0; //mesh->info.real_params[AC_k_pert];
+ // //Wamenumber of the perturbation
+ // const AcReal ampl_pert = xmax/10.0; // xmax/mesh->info.real_params[AC_pert]; //Amplitude of
+ // the perturbation
+ // const AcReal ampl_pert = (AcReal) 0.0;//xmax/20.0; // xmax/mesh->info.real_params[AC_pert];
+ // //Amplitude of the perturbation const AcReal two_pi = (AcReal) 6.28318531;
// const AcReal xorig = mesh->info.real_params[AC_xorig];
// const AcReal zorig = mesh->info.real_params[AC_zorig];
// const AcReal trans = mesh->info.real_params[AC_trans];
-
-
+
// AcReal xx, zz, tanhprof, cosz_wave;
for (int k = 0; k < mz; k++) {
@@ -139,21 +139,19 @@ lnrho_step(AcMesh* mesh)
for (int i = 0; i < mx; i++) {
int idx = i + j * mx + k * mx * my;
// zz = DZ * AcReal(k) - zorig; // Not used
- // cosz_wave = ampl_pert*AcReal(cos(k_pert*((zz/zmax)*two_pi))); // Not used
+ // cosz_wave = ampl_pert*AcReal(cos(k_pert*((zz/zmax)*two_pi))); // Not used
// xx = DX * AcReal(i) - xorig + cosz_wave; //ADD WAVE TODO // Not used
- // tanhprof = AcReal(0.5)*((rho2+rho1) + (rho2-rho1)*AcReal(tanh(xx/trans))); // Not used
- // Commented out the step function initial codition.
- //mesh->vertex_buffer[VTXBUF_LNRHO][idx] = log(tanhprof);
+ // tanhprof = AcReal(0.5)*((rho2+rho1) + (rho2-rho1)*AcReal(tanh(xx/trans))); // Not
+ // used Commented out the step function initial codition.
+ // mesh->vertex_buffer[VTXBUF_LNRHO][idx] = log(tanhprof);
mesh->vertex_buffer[VTXBUF_LNRHO][idx] = lnrho2;
}
}
- }
-
-
+ }
}
// This is the initial condition type for the infalling vedge in the pseudodisk
-// model.
+// model.
void
inflow_vedge(AcMesh* mesh)
{
@@ -170,7 +168,7 @@ inflow_vedge(AcMesh* mesh)
// const double DX = mesh->info.real_params[AC_dsx];
// const double DY = mesh->info.real_params[AC_dsy];
- const double DZ = mesh->info.real_params[AC_dsz];
+ const double DZ = mesh->info.real_params[AC_dsz];
const double AMPL_UU = mesh->info.real_params[AC_ampl_uu];
const double ANGL_UU = mesh->info.real_params[AC_angl_uu];
@@ -184,30 +182,33 @@ inflow_vedge(AcMesh* mesh)
// const AcReal zmax = AcReal(DZ * (nz_max - nz_min));
// const AcReal gaussr = zmax / AcReal(4.0);
- //for (int k = nz_min; k < nz_max; k++) {
+ // for (int k = nz_min; k < nz_max; k++) {
// for (int j = ny_min; j < ny_max; j++) {
// for (int i = nx_min; i < nx_max; i++) {
for (int k = 0; k < mz; k++) {
for (int j = 0; j < my; j++) {
for (int i = 0; i < mx; i++) {
int idx = i + j * mx + k * mx * my;
- zz = DZ * double(k) - zorig;
- //mesh->vertex_buffer[VTXBUF_UUX][idx] = -AMPL_UU*cos(ANGL_UU);
- mesh->vertex_buffer[VTXBUF_UUX][idx] = AcReal(-AMPL_UU*cos(ANGL_UU)*fabs(tanh(zz/trans)));
+ zz = DZ * double(k) - zorig;
+ // mesh->vertex_buffer[VTXBUF_UUX][idx] = -AMPL_UU*cos(ANGL_UU);
+ mesh->vertex_buffer[VTXBUF_UUX][idx] = AcReal(-AMPL_UU * cos(ANGL_UU) *
+ fabs(tanh(zz / trans)));
mesh->vertex_buffer[VTXBUF_UUY][idx] = AcReal(0.0);
- mesh->vertex_buffer[VTXBUF_UUZ][idx] = AcReal(-AMPL_UU*sin(ANGL_UU)*tanh(zz/trans));
+ mesh->vertex_buffer[VTXBUF_UUZ][idx] = AcReal(-AMPL_UU * sin(ANGL_UU) *
+ tanh(zz / trans));
- //Variarion to density
- //AcReal rho = exp(mesh->vertex_buffer[VTXBUF_LNRHO][idx]);
- //NO GAUSSIAN//rho = rho*exp(-(zz/gaussr)*(zz/gaussr));
- //mesh->vertex_buffer[VTXBUF_LNRHO][idx] = log(rho + (range*rho) * (randr() - AcReal(-0.5)));
+ // Variarion to density
+ // AcReal rho = exp(mesh->vertex_buffer[VTXBUF_LNRHO][idx]);
+ // NO GAUSSIAN//rho = rho*exp(-(zz/gaussr)*(zz/gaussr));
+ // mesh->vertex_buffer[VTXBUF_LNRHO][idx] = log(rho + (range*rho) * (randr() -
+ // AcReal(-0.5)));
}
}
}
}
// This is the initial condition type for the infalling vedge in the pseudodisk
-// model.
+// model.
void
inflow_vedge_freefall(AcMesh* mesh)
{
@@ -222,13 +223,13 @@ inflow_vedge_freefall(AcMesh* mesh)
// const int nz_min = mesh->info.int_params[AC_nz_min];
// const int nz_max = mesh->info.int_params[AC_nz_max];
- const double DX = mesh->info.real_params[AC_dsx];
+ const double DX = mesh->info.real_params[AC_dsx];
// const double DY = mesh->info.real_params[AC_dsy];
- const double DZ = mesh->info.real_params[AC_dsz];
+ const double DZ = mesh->info.real_params[AC_dsz];
// const double AMPL_UU = mesh->info.real_params[AC_ampl_uu];
const double ANGL_UU = mesh->info.real_params[AC_angl_uu];
- const double SQ2GM = mesh->info.real_params[AC_sq2GM_star];
+ const double SQ2GM = mesh->info.real_params[AC_sq2GM_star];
// const double GM = mesh->info.real_params[AC_GM_star];
// const double M_star = mesh->info.real_params[AC_M_star];
// const double G_CONST = mesh->info.real_params[AC_G_CONST];
@@ -255,46 +256,51 @@ inflow_vedge_freefall(AcMesh* mesh)
for (int j = 0; j < my; j++) {
for (int i = 0; i < mx; i++) {
int idx = i + j * mx + k * mx * my;
- xx = DX * double(i) - xorig;
- zz = DZ * double(k) - zorig;
+ xx = DX * double(i) - xorig;
+ zz = DZ * double(k) - zorig;
- delx = xx - star_pos_x;
+ delx = xx - star_pos_x;
delz = zz - star_pos_z;
- //TODO: Figure out isthis needed. Now a placeholder.
- //tanhz = fabs(tanh(zz/trans));
+ // TODO: Figure out isthis needed. Now a placeholder.
+ // tanhz = fabs(tanh(zz/trans));
tanhz = 1.0;
-
- RR = sqrt(delx*delx + delz*delz);
- veltot = SQ2GM/sqrt(RR); //Free fall velocity
- //Normal velocity components
- u_x = - veltot*(delx/RR);
- u_z = - veltot*(delz/RR);
+ RR = sqrt(delx * delx + delz * delz);
+ veltot = SQ2GM / sqrt(RR); // Free fall velocity
- //printf("star_pos_z %e, zz %e, delz %e, RR %e\n", star_pos_z, zz, delz, RR);
+ // Normal velocity components
+ u_x = -veltot * (delx / RR);
+ u_z = -veltot * (delz / RR);
- //printf("unit_length = %e, unit_density = %e, unit_velocity = %e,\n M_star = %e, G_CONST = %e, GM = %e, SQ2GM = %e, \n RR = %e, u_x = %e, u_z %e\n",
- // unit_length, unit_density,
+ // printf("star_pos_z %e, zz %e, delz %e, RR %e\n", star_pos_z, zz, delz, RR);
+
+ // printf("unit_length = %e, unit_density = %e, unit_velocity = %e,\n M_star = %e,
+ // G_CONST = %e, GM = %e, SQ2GM = %e, \n RR = %e, u_x = %e, u_z %e\n",
+ // unit_length, unit_density,
// unit_velocity, M_star, G_CONST, GM, SQ2GM, RR, u_x, u_z);
- //printf("%e\n", unit_length*unit_length*unit_length);
+ // printf("%e\n", unit_length*unit_length*unit_length);
-
- //Here including an angel tilt due to pseudodisk
+ // Here including an angel tilt due to pseudodisk
if (delz >= 0.0) {
- mesh->vertex_buffer[VTXBUF_UUX][idx] = AcReal((u_x*cos(ANGL_UU) - u_z*sin(ANGL_UU))*tanhz);
+ mesh->vertex_buffer[VTXBUF_UUX][idx] = AcReal(
+ (u_x * cos(ANGL_UU) - u_z * sin(ANGL_UU)) * tanhz);
mesh->vertex_buffer[VTXBUF_UUY][idx] = AcReal(0.0);
- mesh->vertex_buffer[VTXBUF_UUZ][idx] = AcReal((u_x*sin(ANGL_UU) + u_z*cos(ANGL_UU))*tanhz);
- } else {
- mesh->vertex_buffer[VTXBUF_UUX][idx] = AcReal((u_x*cos(ANGL_UU) + u_z*sin(ANGL_UU))*tanhz);
+ mesh->vertex_buffer[VTXBUF_UUZ][idx] = AcReal(
+ (u_x * sin(ANGL_UU) + u_z * cos(ANGL_UU)) * tanhz);
+ }
+ else {
+ mesh->vertex_buffer[VTXBUF_UUX][idx] = AcReal(
+ (u_x * cos(ANGL_UU) + u_z * sin(ANGL_UU)) * tanhz);
mesh->vertex_buffer[VTXBUF_UUY][idx] = AcReal(0.0);
- mesh->vertex_buffer[VTXBUF_UUZ][idx] = AcReal((-u_x*sin(ANGL_UU) + u_z*cos(ANGL_UU))*tanhz);
+ mesh->vertex_buffer[VTXBUF_UUZ][idx] = AcReal(
+ (-u_x * sin(ANGL_UU) + u_z * cos(ANGL_UU)) * tanhz);
}
}
}
}
}
-// Only x-direction free fall
+// Only x-direction free fall
void
inflow_freefall_x(AcMesh* mesh)
{
@@ -302,7 +308,7 @@ inflow_freefall_x(AcMesh* mesh)
const int my = mesh->info.int_params[AC_my];
const int mz = mesh->info.int_params[AC_mz];
- const double DX = mesh->info.real_params[AC_dsx];
+ const double DX = mesh->info.real_params[AC_dsx];
const double SQ2GM = mesh->info.real_params[AC_sq2GM_star];
// const double G_CONST = mesh->info.real_params[AC_G_CONST];
@@ -320,37 +326,37 @@ inflow_freefall_x(AcMesh* mesh)
for (int j = 0; j < my; j++) {
for (int i = 0; i < mx; i++) {
int idx = i + j * mx + k * mx * my;
- xx = DX * double(i) - xorig;
+ xx = DX * double(i) - xorig;
delx = xx - star_pos_x;
-
+
RR = fabs(delx);
- veltot = SQ2GM/sqrt(RR); //Free fall velocity
+ veltot = SQ2GM / sqrt(RR); // Free fall velocity
- if (isinf(veltot) == 1) printf("xx %e star_pos_x %e delz %e RR %e veltot %e\n",xx, star_pos_x, delx, RR, veltot);
+ if (isinf(veltot) == 1)
+ printf("xx %e star_pos_x %e delz %e RR %e veltot %e\n", xx, star_pos_x, delx,
+ RR, veltot);
- //Normal velocity components
- // u_x = - veltot; // Not used
+ // Normal velocity components
+ // u_x = - veltot; // Not used
- //Freefall condition
- //mesh->vertex_buffer[VTXBUF_UUX][idx] = u_x;
- //mesh->vertex_buffer[VTXBUF_UUY][idx] = 0.0;
- //mesh->vertex_buffer[VTXBUF_UUZ][idx] = 0.0;
+ // Freefall condition
+ // mesh->vertex_buffer[VTXBUF_UUX][idx] = u_x;
+ // mesh->vertex_buffer[VTXBUF_UUY][idx] = 0.0;
+ // mesh->vertex_buffer[VTXBUF_UUZ][idx] = 0.0;
- //Starting with steady state
- mesh->vertex_buffer[VTXBUF_UUX][idx] = 0.0;
+ // Starting with steady state
+ mesh->vertex_buffer[VTXBUF_UUX][idx] = 0.0;
mesh->vertex_buffer[VTXBUF_UUY][idx] = 0.0;
- mesh->vertex_buffer[VTXBUF_UUZ][idx] = 0.0;
+ mesh->vertex_buffer[VTXBUF_UUZ][idx] = 0.0;
- mesh->vertex_buffer[VTXBUF_LNRHO][idx] = AcReal(ampl_lnrho);
+ mesh->vertex_buffer[VTXBUF_LNRHO][idx] = AcReal(ampl_lnrho);
}
}
}
}
-
-
void
gaussian_radial_explosion(AcMesh* mesh)
{
@@ -368,11 +374,11 @@ gaussian_radial_explosion(AcMesh* mesh)
const int nz_min = mesh->info.int_params[AC_nz_min];
const int nz_max = mesh->info.int_params[AC_nz_max];
- const double DX = mesh->info.real_params[AC_dsx];
- const double DY = mesh->info.real_params[AC_dsy];
- const double DZ = mesh->info.real_params[AC_dsz];
+ const double DX = mesh->info.real_params[AC_dsx];
+ const double DY = mesh->info.real_params[AC_dsy];
+ const double DZ = mesh->info.real_params[AC_dsz];
- const double xorig = double(XORIG) - 0.000001;
+ const double xorig = double(XORIG) - 0.000001;
const double yorig = double(YORIG) - 0.000001;
const double zorig = double(ZORIG) - 0.000001;
@@ -447,8 +453,7 @@ gaussian_radial_explosion(AcMesh* mesh)
//+-
yy_abs = -yy;
phi = 2.0 * M_PI - atan(yy_abs / xx);
- if (phi < (3.0 * M_PI) / 2.0 ||
- phi > (2.0 * M_PI + 1e-6)) {
+ if (phi < (3.0 * M_PI) / 2.0 || phi > (2.0 * M_PI + 1e-6)) {
printf("Explosion PHI WRONG +-: xx = %.3f, yy "
"= %.3f, phi = %.3e/PI, M_PI = %.3e\n",
xx, yy, phi / M_PI, M_PI);
@@ -459,23 +464,20 @@ gaussian_radial_explosion(AcMesh* mesh)
yy_abs = -yy;
xx_abs = -xx;
phi = M_PI + atan(yy_abs / xx_abs);
- if (phi < M_PI ||
- phi > ((3.0 * M_PI) / 2.0 + 1e-6)) {
+ if (phi < M_PI || phi > ((3.0 * M_PI) / 2.0 + 1e-6)) {
printf("Explosion PHI WRONG --: xx = %.3f, yy "
"= %.3f, xx_abs = %.3f, yy_abs = %.3f, "
"phi = %.3e, (3.0*M_PI)/2.0 = %.3e\n",
- xx, yy, xx_abs, yy_abs, phi,
- (3.0 * M_PI) / 2.0);
+ xx, yy, xx_abs, yy_abs, phi, (3.0 * M_PI) / 2.0);
}
}
else {
//++
phi = atan(yy / xx);
if (phi < 0 || phi > M_PI / 2.0) {
- printf(
- "Explosion PHI WRONG --: xx = %.3f, yy = "
- "%.3f, phi = %.3e, (3.0*M_PI)/2.0 = %.3e\n",
- xx, yy, phi, (3.0 * M_PI) / 2.0);
+ printf("Explosion PHI WRONG --: xx = %.3f, yy = "
+ "%.3f, phi = %.3e, (3.0*M_PI)/2.0 = %.3e\n",
+ xx, yy, phi, (3.0 * M_PI) / 2.0);
}
}
}
@@ -502,8 +504,8 @@ gaussian_radial_explosion(AcMesh* mesh)
// the exact centre coordinates uu_radial = AMPL_UU*exp(
// -pow((rr - 4.0*WIDTH_UU),2.0) / (2.0*pow(WIDTH_UU, 2.0))
// ); //TODO: Parametrize the peak location.
- uu_radial = AMPL_UU * exp(-pow((rr - UU_SHELL_R), 2.0) /
- (2.0 * pow(WIDTH_UU, 2.0)));
+ uu_radial = AMPL_UU *
+ exp(-pow((rr - UU_SHELL_R), 2.0) / (2.0 * pow(WIDTH_UU, 2.0)));
}
else {
uu_radial = 0.0; // TODO: There will be a discontinuity in
@@ -537,8 +539,7 @@ acmesh_init_to(const InitType& init_type, AcMesh* mesh)
{
srand(123456789);
-
- const int n = AC_VTXBUF_SIZE(mesh->info);
+ const int n = acVertexBufferSize(mesh->info);
const int mx = mesh->info.int_params[AC_mx];
const int my = mesh->info.int_params[AC_my];
@@ -563,7 +564,7 @@ acmesh_init_to(const InitType& init_type, AcMesh* mesh)
}
case INIT_TYPE_GAUSSIAN_RADIAL_EXPL:
acmesh_clear(mesh);
- //acmesh_init_to(INIT_TYPE_RANDOM, mesh);
+ // acmesh_init_to(INIT_TYPE_RANDOM, mesh);
gaussian_radial_explosion(mesh);
break;
@@ -573,21 +574,22 @@ acmesh_init_to(const InitType& init_type, AcMesh* mesh)
for (int k = 0; k < mz; k++) {
for (int j = 0; j < my; j++) {
for (int i = 0; i < mx; i++) {
- int idx = i + j * mx + k * mx * my;
- mesh->vertex_buffer[VTXBUF_UUX][idx] = 2*AcReal(sin(j * AcReal(M_PI) / mx)) - 1;
+ int idx = i + j * mx + k * mx * my;
+ mesh->vertex_buffer[VTXBUF_UUX][idx] = 2 * AcReal(sin(j * AcReal(M_PI) / mx)) -
+ 1;
}
}
}
break;
- case INIT_TYPE_VEDGE:
+ case INIT_TYPE_VEDGE:
acmesh_clear(mesh);
inflow_vedge_freefall(mesh);
break;
- case INIT_TYPE_VEDGEX:
+ case INIT_TYPE_VEDGEX:
acmesh_clear(mesh);
inflow_freefall_x(mesh);
break;
- case INIT_TYPE_RAYLEIGH_TAYLOR:
+ case INIT_TYPE_RAYLEIGH_TAYLOR:
acmesh_clear(mesh);
inflow_freefall_x(mesh);
lnrho_step(mesh);
@@ -627,9 +629,15 @@ acmesh_init_to(const InitType& init_type, AcMesh* mesh)
const AcReal ky_uu = 8.;
const AcReal kz_uu = 8.;
- mesh->vertex_buffer[VTXBUF_UUX][idx] = ampl_uu * (ABC_A * (AcReal)sin(kz_uu * zz) + ABC_C * (AcReal)cos(ky_uu * yy));
- mesh->vertex_buffer[VTXBUF_UUY][idx] = ampl_uu * (ABC_B * (AcReal)sin(kx_uu * xx) + ABC_A * (AcReal)cos(kz_uu * zz));
- mesh->vertex_buffer[VTXBUF_UUZ][idx] = ampl_uu * (ABC_C * (AcReal)sin(ky_uu * yy) + ABC_B * (AcReal)cos(kx_uu * xx));
+ mesh->vertex_buffer[VTXBUF_UUX][idx] = ampl_uu *
+ (ABC_A * (AcReal)sin(kz_uu * zz) +
+ ABC_C * (AcReal)cos(ky_uu * yy));
+ mesh->vertex_buffer[VTXBUF_UUY][idx] = ampl_uu *
+ (ABC_B * (AcReal)sin(kx_uu * xx) +
+ ABC_A * (AcReal)cos(kz_uu * zz));
+ mesh->vertex_buffer[VTXBUF_UUZ][idx] = ampl_uu *
+ (ABC_C * (AcReal)sin(ky_uu * yy) +
+ ABC_B * (AcReal)cos(kx_uu * xx));
}
}
}
@@ -637,20 +645,23 @@ acmesh_init_to(const InitType& init_type, AcMesh* mesh)
}
case INIT_TYPE_RAYLEIGH_BENARD: {
acmesh_init_to(INIT_TYPE_RANDOM, mesh);
- #if LTEMPERATURE
+#if LTEMPERATURE
vertex_buffer_set(VTXBUF_LNRHO, 1, mesh);
const AcReal range = AcReal(0.9);
for (int k = nz_min; k < nz_max; k++) {
for (int j = ny_min; j < ny_max; j++) {
for (int i = nx_min; i < nx_max; i++) {
- const int idx = i + j * mx + k * mx * my;
- mesh->vertex_buffer[VTXBUF_TEMPERATURE][idx] = (range * (k - nz_min)) / mesh->info.int_params[AC_nz] + 0.1;
+ const int idx = i + j * mx + k * mx * my;
+ mesh->vertex_buffer[VTXBUF_TEMPERATURE][idx] = (range * (k - nz_min)) /
+ mesh->info
+ .int_params[AC_nz] +
+ 0.1;
}
}
}
- #else
+#else
WARNING("INIT_TYPE_RAYLEIGH_BERNARD called even though VTXBUF_TEMPERATURE is not used");
- #endif
+#endif
break;
}
default:
@@ -688,14 +699,13 @@ acmesh_destroy(AcMesh* mesh)
free(mesh);
}
-
ModelMesh*
modelmesh_create(const AcMeshInfo& mesh_info)
{
ModelMesh* mesh = (ModelMesh*)malloc(sizeof(*mesh));
- mesh->info = mesh_info;
+ mesh->info = mesh_info;
- const size_t bytes = AC_VTXBUF_SIZE(mesh->info) * sizeof(mesh->vertex_buffer[0][0]);
+ const size_t bytes = acVertexBufferSize(mesh->info) * sizeof(mesh->vertex_buffer[0][0]);
for (int i = 0; i < NUM_VTXBUF_HANDLES; ++i) {
mesh->vertex_buffer[VertexBufferHandle(i)] = (ModelScalar*)malloc(bytes);
ERRCHK(mesh->vertex_buffer[VertexBufferHandle(i)] != NULL);
@@ -721,7 +731,7 @@ acmesh_to_modelmesh(const AcMesh& acmesh, ModelMesh* modelmesh)
memcpy(&modelmesh->info, &acmesh.info, sizeof(acmesh.info));
for (int i = 0; i < NUM_VTXBUF_HANDLES; ++i)
- for (size_t j = 0; j < AC_VTXBUF_SIZE(acmesh.info); ++j)
+ for (size_t j = 0; j < acVertexBufferSize(acmesh.info); ++j)
modelmesh->vertex_buffer[i][j] = (ModelScalar)acmesh.vertex_buffer[i][j];
}
@@ -732,6 +742,6 @@ modelmesh_to_acmesh(const ModelMesh& modelmesh, AcMesh* acmesh)
memcpy(&acmesh->info, &modelmesh.info, sizeof(modelmesh.info));
for (int i = 0; i < NUM_VTXBUF_HANDLES; ++i)
- for (size_t j = 0; j < AC_VTXBUF_SIZE(modelmesh.info); ++j)
+ for (size_t j = 0; j < acVertexBufferSize(modelmesh.info); ++j)
acmesh->vertex_buffer[i][j] = (AcReal)modelmesh.vertex_buffer[i][j];
}
diff --git a/src/standalone/model/host_memory.h b/src/standalone/model/host_memory.h
index 105396c..e2958cd 100644
--- a/src/standalone/model/host_memory.h
+++ b/src/standalone/model/host_memory.h
@@ -40,7 +40,9 @@
FUNC(INIT_TYPE_RAYLEIGH_BENARD)
// clang-format on
+#define AC_GEN_ID(X) X
typedef enum { AC_FOR_INIT_TYPES(AC_GEN_ID), NUM_INIT_TYPES } InitType;
+#undef AC_GEN_ID
extern const char* init_type_names[]; // Defined in host_memory.cc
diff --git a/src/standalone/model/model_boundconds.cc b/src/standalone/model/model_boundconds.cc
index 5d107ae..188b97e 100644
--- a/src/standalone/model/model_boundconds.cc
+++ b/src/standalone/model/model_boundconds.cc
@@ -86,10 +86,10 @@ boundconds(const AcMeshInfo& mesh_info, ModelMesh* mesh)
j_src += ny_min;
k_src += nz_min;
- const size_t src_idx = AC_VTXBUF_IDX(i_src, j_src, k_src, mesh_info);
- const size_t dst_idx = AC_VTXBUF_IDX(i_dst, j_dst, k_dst, mesh_info);
- ERRCHK(src_idx < AC_VTXBUF_SIZE(mesh_info));
- ERRCHK(dst_idx < AC_VTXBUF_SIZE(mesh_info));
+ const size_t src_idx = acVertexBufferIdx(i_src, j_src, k_src, mesh_info);
+ const size_t dst_idx = acVertexBufferIdx(i_dst, j_dst, k_dst, mesh_info);
+ ERRCHK(src_idx < acVertexBufferSize(mesh_info));
+ ERRCHK(dst_idx < acVertexBufferSize(mesh_info));
mesh->vertex_buffer[w][dst_idx] = mesh->vertex_buffer[w][src_idx];
}
}
diff --git a/src/standalone/model/model_diff.h b/src/standalone/model/model_diff.h
index 186c984..20678bf 100644
--- a/src/standalone/model/model_diff.h
+++ b/src/standalone/model/model_diff.h
@@ -41,30 +41,30 @@ der_scal(const int& i, const int& j, const int& k, const AcMeshInfo& mesh_info,
switch (axis) {
case AXIS_X:
- f0 = scal[AC_VTXBUF_IDX(i - 3, j, k, mesh_info)];
- f1 = scal[AC_VTXBUF_IDX(i - 2, j, k, mesh_info)];
- f2 = scal[AC_VTXBUF_IDX(i - 1, j, k, mesh_info)];
- f4 = scal[AC_VTXBUF_IDX(i + 1, j, k, mesh_info)];
- f5 = scal[AC_VTXBUF_IDX(i + 2, j, k, mesh_info)];
- f6 = scal[AC_VTXBUF_IDX(i + 3, j, k, mesh_info)];
+ f0 = scal[acVertexBufferIdx(i - 3, j, k, mesh_info)];
+ f1 = scal[acVertexBufferIdx(i - 2, j, k, mesh_info)];
+ f2 = scal[acVertexBufferIdx(i - 1, j, k, mesh_info)];
+ f4 = scal[acVertexBufferIdx(i + 1, j, k, mesh_info)];
+ f5 = scal[acVertexBufferIdx(i + 2, j, k, mesh_info)];
+ f6 = scal[acVertexBufferIdx(i + 3, j, k, mesh_info)];
ds = mesh_info.real_params[AC_dsx];
break;
case AXIS_Y:
- f0 = scal[AC_VTXBUF_IDX(i, j - 3, k, mesh_info)];
- f1 = scal[AC_VTXBUF_IDX(i, j - 2, k, mesh_info)];
- f2 = scal[AC_VTXBUF_IDX(i, j - 1, k, mesh_info)];
- f4 = scal[AC_VTXBUF_IDX(i, j + 1, k, mesh_info)];
- f5 = scal[AC_VTXBUF_IDX(i, j + 2, k, mesh_info)];
- f6 = scal[AC_VTXBUF_IDX(i, j + 3, k, mesh_info)];
+ f0 = scal[acVertexBufferIdx(i, j - 3, k, mesh_info)];
+ f1 = scal[acVertexBufferIdx(i, j - 2, k, mesh_info)];
+ f2 = scal[acVertexBufferIdx(i, j - 1, k, mesh_info)];
+ f4 = scal[acVertexBufferIdx(i, j + 1, k, mesh_info)];
+ f5 = scal[acVertexBufferIdx(i, j + 2, k, mesh_info)];
+ f6 = scal[acVertexBufferIdx(i, j + 3, k, mesh_info)];
ds = mesh_info.real_params[AC_dsy];
break;
case AXIS_Z:
- f0 = scal[AC_VTXBUF_IDX(i, j, k - 3, mesh_info)];
- f1 = scal[AC_VTXBUF_IDX(i, j, k - 2, mesh_info)];
- f2 = scal[AC_VTXBUF_IDX(i, j, k - 1, mesh_info)];
- f4 = scal[AC_VTXBUF_IDX(i, j, k + 1, mesh_info)];
- f5 = scal[AC_VTXBUF_IDX(i, j, k + 2, mesh_info)];
- f6 = scal[AC_VTXBUF_IDX(i, j, k + 3, mesh_info)];
+ f0 = scal[acVertexBufferIdx(i, j, k - 3, mesh_info)];
+ f1 = scal[acVertexBufferIdx(i, j, k - 2, mesh_info)];
+ f2 = scal[acVertexBufferIdx(i, j, k - 1, mesh_info)];
+ f4 = scal[acVertexBufferIdx(i, j, k + 1, mesh_info)];
+ f5 = scal[acVertexBufferIdx(i, j, k + 2, mesh_info)];
+ f6 = scal[acVertexBufferIdx(i, j, k + 3, mesh_info)];
ds = mesh_info.real_params[AC_dsz];
break;
default:
@@ -82,34 +82,34 @@ der2_scal(const int& i, const int& j, const int& k, const AcMeshInfo& mesh_info,
MODEL_REAL f0, f1, f2, f3, f4, f5, f6;
MODEL_REAL ds;
- f3 = scal[AC_VTXBUF_IDX(i, j, k, mesh_info)];
+ f3 = scal[acVertexBufferIdx(i, j, k, mesh_info)];
switch (axis) {
case AXIS_X:
- f0 = scal[AC_VTXBUF_IDX(i - 3, j, k, mesh_info)];
- f1 = scal[AC_VTXBUF_IDX(i - 2, j, k, mesh_info)];
- f2 = scal[AC_VTXBUF_IDX(i - 1, j, k, mesh_info)];
- f4 = scal[AC_VTXBUF_IDX(i + 1, j, k, mesh_info)];
- f5 = scal[AC_VTXBUF_IDX(i + 2, j, k, mesh_info)];
- f6 = scal[AC_VTXBUF_IDX(i + 3, j, k, mesh_info)];
+ f0 = scal[acVertexBufferIdx(i - 3, j, k, mesh_info)];
+ f1 = scal[acVertexBufferIdx(i - 2, j, k, mesh_info)];
+ f2 = scal[acVertexBufferIdx(i - 1, j, k, mesh_info)];
+ f4 = scal[acVertexBufferIdx(i + 1, j, k, mesh_info)];
+ f5 = scal[acVertexBufferIdx(i + 2, j, k, mesh_info)];
+ f6 = scal[acVertexBufferIdx(i + 3, j, k, mesh_info)];
ds = mesh_info.real_params[AC_dsx];
break;
case AXIS_Y:
- f0 = scal[AC_VTXBUF_IDX(i, j - 3, k, mesh_info)];
- f1 = scal[AC_VTXBUF_IDX(i, j - 2, k, mesh_info)];
- f2 = scal[AC_VTXBUF_IDX(i, j - 1, k, mesh_info)];
- f4 = scal[AC_VTXBUF_IDX(i, j + 1, k, mesh_info)];
- f5 = scal[AC_VTXBUF_IDX(i, j + 2, k, mesh_info)];
- f6 = scal[AC_VTXBUF_IDX(i, j + 3, k, mesh_info)];
+ f0 = scal[acVertexBufferIdx(i, j - 3, k, mesh_info)];
+ f1 = scal[acVertexBufferIdx(i, j - 2, k, mesh_info)];
+ f2 = scal[acVertexBufferIdx(i, j - 1, k, mesh_info)];
+ f4 = scal[acVertexBufferIdx(i, j + 1, k, mesh_info)];
+ f5 = scal[acVertexBufferIdx(i, j + 2, k, mesh_info)];
+ f6 = scal[acVertexBufferIdx(i, j + 3, k, mesh_info)];
ds = mesh_info.real_params[AC_dsy];
break;
case AXIS_Z:
- f0 = scal[AC_VTXBUF_IDX(i, j, k - 3, mesh_info)];
- f1 = scal[AC_VTXBUF_IDX(i, j, k - 2, mesh_info)];
- f2 = scal[AC_VTXBUF_IDX(i, j, k - 1, mesh_info)];
- f4 = scal[AC_VTXBUF_IDX(i, j, k + 1, mesh_info)];
- f5 = scal[AC_VTXBUF_IDX(i, j, k + 2, mesh_info)];
- f6 = scal[AC_VTXBUF_IDX(i, j, k + 3, mesh_info)];
+ f0 = scal[acVertexBufferIdx(i, j, k - 3, mesh_info)];
+ f1 = scal[acVertexBufferIdx(i, j, k - 2, mesh_info)];
+ f2 = scal[acVertexBufferIdx(i, j, k - 1, mesh_info)];
+ f4 = scal[acVertexBufferIdx(i, j, k + 1, mesh_info)];
+ f5 = scal[acVertexBufferIdx(i, j, k + 2, mesh_info)];
+ f6 = scal[acVertexBufferIdx(i, j, k + 3, mesh_info)];
ds = mesh_info.real_params[AC_dsz];
break;
default:
@@ -163,7 +163,7 @@ vec_dot_nabla_scal(const int& i, const int& j, const int& k,
const AcMeshInfo& mesh_info, const MODEL_REAL* vec_x,
const MODEL_REAL* vec_y, const MODEL_REAL* vec_z, const MODEL_REAL* scal)
{
- const int idx = AC_VTXBUF_IDX(i, j, k, mesh_info);
+ const int idx = acVertexBufferIdx(i, j, k, mesh_info);
MODEL_REAL ddx_scal, ddy_scal, ddz_scal;
grad(i, j, k, mesh_info, scal, &ddx_scal, &ddy_scal, &ddz_scal);
return vec_x[idx] * ddx_scal + vec_y[idx] * ddy_scal +
@@ -196,56 +196,56 @@ dernm_scal(const int& i, const int& j, const int& k,
switch (dernm) {
case DERNM_XY:
fac = MODEL_REAL(1. / 720.) * (MODEL_REAL(1.) / dsx) * (MODEL_REAL(1.) / dsy);
- f_p1_p1 = scal[AC_VTXBUF_IDX(i + 1, j + 1, k, mesh_info)];
- f_m1_p1 = scal[AC_VTXBUF_IDX(i - 1, j + 1, k, mesh_info)];
- f_m1_m1 = scal[AC_VTXBUF_IDX(i - 1, j - 1, k, mesh_info)];
- f_p1_m1 = scal[AC_VTXBUF_IDX(i + 1, j - 1, k, mesh_info)];
+ f_p1_p1 = scal[acVertexBufferIdx(i + 1, j + 1, k, mesh_info)];
+ f_m1_p1 = scal[acVertexBufferIdx(i - 1, j + 1, k, mesh_info)];
+ f_m1_m1 = scal[acVertexBufferIdx(i - 1, j - 1, k, mesh_info)];
+ f_p1_m1 = scal[acVertexBufferIdx(i + 1, j - 1, k, mesh_info)];
- f_p2_p2 = scal[AC_VTXBUF_IDX(i + 2, j + 2, k, mesh_info)];
- f_m2_p2 = scal[AC_VTXBUF_IDX(i - 2, j + 2, k, mesh_info)];
- f_m2_m2 = scal[AC_VTXBUF_IDX(i - 2, j - 2, k, mesh_info)];
- f_p2_m2 = scal[AC_VTXBUF_IDX(i + 2, j - 2, k, mesh_info)];
+ f_p2_p2 = scal[acVertexBufferIdx(i + 2, j + 2, k, mesh_info)];
+ f_m2_p2 = scal[acVertexBufferIdx(i - 2, j + 2, k, mesh_info)];
+ f_m2_m2 = scal[acVertexBufferIdx(i - 2, j - 2, k, mesh_info)];
+ f_p2_m2 = scal[acVertexBufferIdx(i + 2, j - 2, k, mesh_info)];
- f_p3_p3 = scal[AC_VTXBUF_IDX(i + 3, j + 3, k, mesh_info)];
- f_m3_p3 = scal[AC_VTXBUF_IDX(i - 3, j + 3, k, mesh_info)];
- f_m3_m3 = scal[AC_VTXBUF_IDX(i - 3, j - 3, k, mesh_info)];
- f_p3_m3 = scal[AC_VTXBUF_IDX(i + 3, j - 3, k, mesh_info)];
+ f_p3_p3 = scal[acVertexBufferIdx(i + 3, j + 3, k, mesh_info)];
+ f_m3_p3 = scal[acVertexBufferIdx(i - 3, j + 3, k, mesh_info)];
+ f_m3_m3 = scal[acVertexBufferIdx(i - 3, j - 3, k, mesh_info)];
+ f_p3_m3 = scal[acVertexBufferIdx(i + 3, j - 3, k, mesh_info)];
break;
case DERNM_YZ:
// NOTE this is a bit different from the old one, second is j+1k-1
// instead of j-1,k+1
fac = MODEL_REAL(1. / 720.) * (MODEL_REAL(1.) / dsy) * (MODEL_REAL(1.) / dsz);
- f_p1_p1 = scal[AC_VTXBUF_IDX(i, j + 1, k + 1, mesh_info)];
- f_m1_p1 = scal[AC_VTXBUF_IDX(i, j - 1, k + 1, mesh_info)];
- f_m1_m1 = scal[AC_VTXBUF_IDX(i, j - 1, k - 1, mesh_info)];
- f_p1_m1 = scal[AC_VTXBUF_IDX(i, j + 1, k - 1, mesh_info)];
+ f_p1_p1 = scal[acVertexBufferIdx(i, j + 1, k + 1, mesh_info)];
+ f_m1_p1 = scal[acVertexBufferIdx(i, j - 1, k + 1, mesh_info)];
+ f_m1_m1 = scal[acVertexBufferIdx(i, j - 1, k - 1, mesh_info)];
+ f_p1_m1 = scal[acVertexBufferIdx(i, j + 1, k - 1, mesh_info)];
- f_p2_p2 = scal[AC_VTXBUF_IDX(i, j + 2, k + 2, mesh_info)];
- f_m2_p2 = scal[AC_VTXBUF_IDX(i, j - 2, k + 2, mesh_info)];
- f_m2_m2 = scal[AC_VTXBUF_IDX(i, j - 2, k - 2, mesh_info)];
- f_p2_m2 = scal[AC_VTXBUF_IDX(i, j + 2, k - 2, mesh_info)];
+ f_p2_p2 = scal[acVertexBufferIdx(i, j + 2, k + 2, mesh_info)];
+ f_m2_p2 = scal[acVertexBufferIdx(i, j - 2, k + 2, mesh_info)];
+ f_m2_m2 = scal[acVertexBufferIdx(i, j - 2, k - 2, mesh_info)];
+ f_p2_m2 = scal[acVertexBufferIdx(i, j + 2, k - 2, mesh_info)];
- f_p3_p3 = scal[AC_VTXBUF_IDX(i, j + 3, k + 3, mesh_info)];
- f_m3_p3 = scal[AC_VTXBUF_IDX(i, j - 3, k + 3, mesh_info)];
- f_m3_m3 = scal[AC_VTXBUF_IDX(i, j - 3, k - 3, mesh_info)];
- f_p3_m3 = scal[AC_VTXBUF_IDX(i, j + 3, k - 3, mesh_info)];
+ f_p3_p3 = scal[acVertexBufferIdx(i, j + 3, k + 3, mesh_info)];
+ f_m3_p3 = scal[acVertexBufferIdx(i, j - 3, k + 3, mesh_info)];
+ f_m3_m3 = scal[acVertexBufferIdx(i, j - 3, k - 3, mesh_info)];
+ f_p3_m3 = scal[acVertexBufferIdx(i, j + 3, k - 3, mesh_info)];
break;
case DERNM_XZ:
fac = MODEL_REAL(1. / 720.) * (MODEL_REAL(1.) / dsx) * (MODEL_REAL(1.) / dsz);
- f_p1_p1 = scal[AC_VTXBUF_IDX(i + 1, j, k + 1, mesh_info)];
- f_m1_p1 = scal[AC_VTXBUF_IDX(i - 1, j, k + 1, mesh_info)];
- f_m1_m1 = scal[AC_VTXBUF_IDX(i - 1, j, k - 1, mesh_info)];
- f_p1_m1 = scal[AC_VTXBUF_IDX(i + 1, j, k - 1, mesh_info)];
+ f_p1_p1 = scal[acVertexBufferIdx(i + 1, j, k + 1, mesh_info)];
+ f_m1_p1 = scal[acVertexBufferIdx(i - 1, j, k + 1, mesh_info)];
+ f_m1_m1 = scal[acVertexBufferIdx(i - 1, j, k - 1, mesh_info)];
+ f_p1_m1 = scal[acVertexBufferIdx(i + 1, j, k - 1, mesh_info)];
- f_p2_p2 = scal[AC_VTXBUF_IDX(i + 2, j, k + 2, mesh_info)];
- f_m2_p2 = scal[AC_VTXBUF_IDX(i - 2, j, k + 2, mesh_info)];
- f_m2_m2 = scal[AC_VTXBUF_IDX(i - 2, j, k - 2, mesh_info)];
- f_p2_m2 = scal[AC_VTXBUF_IDX(i + 2, j, k - 2, mesh_info)];
+ f_p2_p2 = scal[acVertexBufferIdx(i + 2, j, k + 2, mesh_info)];
+ f_m2_p2 = scal[acVertexBufferIdx(i - 2, j, k + 2, mesh_info)];
+ f_m2_m2 = scal[acVertexBufferIdx(i - 2, j, k - 2, mesh_info)];
+ f_p2_m2 = scal[acVertexBufferIdx(i + 2, j, k - 2, mesh_info)];
- f_p3_p3 = scal[AC_VTXBUF_IDX(i + 3, j, k + 3, mesh_info)];
- f_m3_p3 = scal[AC_VTXBUF_IDX(i - 3, j, k + 3, mesh_info)];
- f_m3_m3 = scal[AC_VTXBUF_IDX(i - 3, j, k - 3, mesh_info)];
- f_p3_m3 = scal[AC_VTXBUF_IDX(i + 3, j, k - 3, mesh_info)];
+ f_p3_p3 = scal[acVertexBufferIdx(i + 3, j, k + 3, mesh_info)];
+ f_m3_p3 = scal[acVertexBufferIdx(i - 3, j, k + 3, mesh_info)];
+ f_m3_m3 = scal[acVertexBufferIdx(i - 3, j, k - 3, mesh_info)];
+ f_p3_m3 = scal[acVertexBufferIdx(i + 3, j, k - 3, mesh_info)];
break;
default:
ERROR("Invalid dernm type");
diff --git a/src/standalone/model/model_reduce.cc b/src/standalone/model/model_reduce.cc
index dde5403..6d48c4b 100644
--- a/src/standalone/model/model_reduce.cc
+++ b/src/standalone/model/model_reduce.cc
@@ -99,7 +99,7 @@ model_reduce_scal(const ModelMesh& mesh, const ReductionType& rtype,
ERROR("Unrecognized RTYPE");
}
- const int initial_idx = AC_VTXBUF_IDX(
+ const int initial_idx = acVertexBufferIdx(
mesh.info.int_params[AC_nx_min], mesh.info.int_params[AC_ny_min],
mesh.info.int_params[AC_nz_min], mesh.info);
@@ -115,7 +115,7 @@ model_reduce_scal(const ModelMesh& mesh, const ReductionType& rtype,
j < mesh.info.int_params[AC_ny_max]; ++j) {
for (int i = mesh.info.int_params[AC_nx_min];
i < mesh.info.int_params[AC_nx_max]; ++i) {
- const int idx = AC_VTXBUF_IDX(i, j, k, mesh.info);
+ const int idx = acVertexBufferIdx(i, j, k, mesh.info);
const ModelScalar curr_val = reduce_initial(
mesh.vertex_buffer[a][idx]);
res = reduce(res, curr_val);
@@ -166,7 +166,7 @@ model_reduce_vec(const ModelMesh& mesh, const ReductionType& rtype,
ERROR("Unrecognized RTYPE");
}
- const int initial_idx = AC_VTXBUF_IDX(
+ const int initial_idx = acVertexBufferIdx(
mesh.info.int_params[AC_nx_min], mesh.info.int_params[AC_ny_min],
mesh.info.int_params[AC_nz_min], mesh.info);
@@ -184,7 +184,7 @@ model_reduce_vec(const ModelMesh& mesh, const ReductionType& rtype,
j < mesh.info.int_params[AC_ny_max]; j++) {
for (int i = mesh.info.int_params[AC_nx_min];
i < mesh.info.int_params[AC_nx_max]; i++) {
- const int idx = AC_VTXBUF_IDX(i, j, k, mesh.info);
+ const int idx = acVertexBufferIdx(i, j, k, mesh.info);
const ModelScalar curr_val = reduce_initial(
mesh.vertex_buffer[a][idx], mesh.vertex_buffer[b][idx],
mesh.vertex_buffer[c][idx]);
diff --git a/src/standalone/model/model_rk3.cc b/src/standalone/model/model_rk3.cc
index 9353980..3907f93 100644
--- a/src/standalone/model/model_rk3.cc
+++ b/src/standalone/model/model_rk3.cc
@@ -68,7 +68,7 @@ get(const AcRealParam param)
static inline int
IDX(const int i, const int j, const int k)
{
- return AC_VTXBUF_IDX(i, j, k, (*mesh_info));
+ return acVertexBufferIdx(i, j, k, (*mesh_info));
}
/*
@@ -749,7 +749,7 @@ static void
solve_alpha_step(const int step_number, const ModelScalar dt, const int i, const int j, const int k,
ModelMesh& in, ModelMesh* out)
{
- const int idx = AC_VTXBUF_IDX(i, j, k, in.info);
+ const int idx = acVertexBufferIdx(i, j, k, in.info);
const ModelScalarData lnrho = read_data(i, j, k, in.vertex_buffer, VTXBUF_LNRHO);
const ModelVectorData uu = read_data(i, j, k, in.vertex_buffer,
@@ -797,7 +797,7 @@ static void
solve_beta_step(const int step_number, const ModelScalar dt, const int i, const int j, const int k,
const ModelMesh& in, ModelMesh* out)
{
- const int idx = AC_VTXBUF_IDX(i, j, k, in.info);
+ const int idx = acVertexBufferIdx(i, j, k, in.info);
// Williamson (1980) NOTE: older version of astaroth used inhomogenous
const ModelScalar beta[] = {ModelScalar(1. / 3.), ModelScalar(15. / 16.),
diff --git a/src/standalone/renderer.cc b/src/standalone/renderer.cc
index 52aed16..1522cc5 100644
--- a/src/standalone/renderer.cc
+++ b/src/standalone/renderer.cc
@@ -162,7 +162,7 @@ draw_vertex_buffer(const AcMesh& mesh, const VertexBufferHandle& vertex_buffer,
for (int i = 0; i < mesh.info.int_params[AC_mx]; ++i) {
ERRCHK(i < datasurface_width && j < datasurface_height);
- const int idx = AC_VTXBUF_IDX(i, j, k, mesh.info);
+ const int idx = acVertexBufferIdx(i, j, k, mesh.info);
const uint8_t shade = (uint8_t)(
255.f * (fabsf(float(mesh.vertex_buffer[vertex_buffer][idx]) - mid)) / range);
uint8_t color[4] = {0, 0, 0, 255};
@@ -219,7 +219,7 @@ draw_vertex_buffer_vec(const AcMesh& mesh, const VertexBufferHandle& vertex_buff
for (int i = 0; i < mesh.info.int_params[AC_mx]; ++i) {
ERRCHK(i < datasurface_width && j < datasurface_height);
- const int idx = AC_VTXBUF_IDX(i, j, k, mesh.info);
+ const int idx = acVertexBufferIdx(i, j, k, mesh.info);
const uint8_t r = (uint8_t)(
255.f * (fabsf(float(mesh.vertex_buffer[vertex_buffer_a][idx]) - mid)) / range);
const uint8_t g = (uint8_t)(
diff --git a/src/standalone/simulation.cc b/src/standalone/simulation.cc
index af68b2d..c6842dc 100644
--- a/src/standalone/simulation.cc
+++ b/src/standalone/simulation.cc
@@ -100,7 +100,7 @@ save_mesh(const AcMesh& save_mesh, const int step, const AcReal t_step)
FILE* save_ptr;
for (int w = 0; w < NUM_VTXBUF_HANDLES; ++w) {
- const size_t n = AC_VTXBUF_SIZE(save_mesh.info);
+ const size_t n = acVertexBufferSize(save_mesh.info);
const char* buffername = vtxbuf_names[w];
char cstep[11];