Merge branch 'master' into sink_20190723
This commit is contained in:
@@ -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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
#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) \
|
||||
#if LENTROPY
|
||||
#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_HYDRO_VTXBUF_HANDLES(FUNC)
|
||||
#define AC_FOR_VTXBUF_HANDLES(FUNC) \
|
||||
FUNC(VTXBUF_LNRHO),
|
||||
#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),
|
||||
#else
|
||||
#define AC_FOR_ENTROPY_VTXBUF_HANDLES(FUNC)
|
||||
#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
|
||||
|
@@ -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])
|
||||
|
||||
|
@@ -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)]));
|
||||
|
||||
```
|
||||
|
||||
|
@@ -16,209 +16,13 @@
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with Astaroth. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
/**
|
||||
* @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 <float.h> // FLT_EPSILON, etc
|
||||
#include <stdlib.h> // size_t
|
||||
#include <vector_types.h> // 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
|
||||
|
226
include/astaroth_defines.h
Normal file
226
include/astaroth_defines.h
Normal file
@@ -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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
#pragma once
|
||||
|
||||
#ifdef __cplusplus
|
||||
extern "C" {
|
||||
#endif
|
||||
|
||||
#include <float.h> // FLT_EPSILON, etc
|
||||
#include <stdlib.h> // size_t
|
||||
#include <vector_types.h> // 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
|
@@ -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)
|
||||
#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* realparam_names[] = {AC_FOR_REAL_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;
|
||||
}
|
||||
|
@@ -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,
|
||||
|
@@ -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]);
|
||||
|
||||
|
@@ -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));
|
||||
}
|
||||
|
||||
|
@@ -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])
|
||||
@@ -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,7 +95,6 @@ randr(void)
|
||||
return AcReal(rand()) / AcReal(RAND_MAX);
|
||||
}
|
||||
|
||||
|
||||
void
|
||||
lnrho_step(AcMesh* mesh)
|
||||
{
|
||||
@@ -117,21 +116,22 @@ lnrho_step(AcMesh* mesh)
|
||||
// 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 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++) {
|
||||
@@ -141,15 +141,13 @@ lnrho_step(AcMesh* mesh)
|
||||
// zz = DZ * AcReal(k) - zorig; // 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
|
||||
@@ -184,7 +182,7 @@ 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++) {
|
||||
@@ -192,15 +190,18 @@ inflow_vedge(AcMesh* mesh)
|
||||
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)));
|
||||
// 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)));
|
||||
}
|
||||
}
|
||||
}
|
||||
@@ -260,34 +261,39 @@ inflow_vedge_freefall(AcMesh* mesh)
|
||||
|
||||
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
|
||||
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);
|
||||
// Normal velocity components
|
||||
u_x = -veltot * (delx / RR);
|
||||
u_z = -veltot * (delz / RR);
|
||||
|
||||
//printf("star_pos_z %e, zz %e, delz %e, RR %e\n", star_pos_z, zz, delz, RR);
|
||||
// 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",
|
||||
// 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);
|
||||
}
|
||||
}
|
||||
}
|
||||
@@ -326,19 +332,21 @@ inflow_freefall_x(AcMesh* mesh)
|
||||
|
||||
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
|
||||
// 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
|
||||
// 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;
|
||||
@@ -349,8 +357,6 @@ inflow_freefall_x(AcMesh* mesh)
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
void
|
||||
gaussian_radial_explosion(AcMesh* mesh)
|
||||
{
|
||||
@@ -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,21 +464,18 @@ 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 = "
|
||||
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;
|
||||
@@ -574,7 +575,8 @@ acmesh_init_to(const InitType& init_type, AcMesh* mesh)
|
||||
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;
|
||||
mesh->vertex_buffer[VTXBUF_UUX][idx] = 2 * AcReal(sin(j * AcReal(M_PI) / mx)) -
|
||||
1;
|
||||
}
|
||||
}
|
||||
}
|
||||
@@ -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;
|
||||
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;
|
||||
|
||||
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];
|
||||
}
|
||||
|
@@ -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
|
||||
|
||||
|
@@ -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];
|
||||
}
|
||||
}
|
||||
|
@@ -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");
|
||||
|
@@ -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]);
|
||||
|
@@ -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.),
|
||||
|
@@ -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)(
|
||||
|
@@ -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];
|
||||
|
Reference in New Issue
Block a user