diff options
Diffstat (limited to 'src/qms.h')
| -rw-r--r-- | src/qms.h | 301 |
1 files changed, 301 insertions, 0 deletions
diff --git a/src/qms.h b/src/qms.h new file mode 100644 index 0000000..f88c0a8 --- /dev/null +++ b/src/qms.h @@ -0,0 +1,301 @@ + +#include <inttypes.h> + +#include <cl.h> + +#include "cctk.h" + +#define POLAR (1) + +#define CCZ4 (0) + +#define SQR(x) ((x) * (x)) +#define SGN(x) ((x) >= 0.0 ? 1.0 : -1.0) +#define MAX(x, y) ((x) > (y) ? (x) : (y)) +#define MIN(x, y) ((x) > (y) ? (y) : (x)) +#define ARRAY_ELEMS(arr) (sizeof(arr) / sizeof(*arr)) + +/* + * small number to avoid r=0 singularities + */ +#define EPS 1E-08 + +#define SCALE_FACTOR scale_factor + +/* indices (in our code, not cactus structs) of the grid functions which we'll need to + * interpolate on the pseudospectral grid */ +enum MetricVars { + GTXX = 0, + GTYY, + GTZZ, + GTXY, + GTXZ, + GTYZ, + PHI, + ATXX, + ATYY, + ATZZ, + ATXY, + ATXZ, + ATYZ, + K, + XTX, + XTY, + XTZ, + BETAX, + BETAY, + BETAZ, + ALPHA, + KDOT_XX, + KDOT_YY, + KDOT_ZZ, + KDOT_XY, + KDOT_XZ, + KDOT_YZ, + XTDOT_X, + XTDOT_Y, + XTDOT_Z, + PHIDOT, + NB_METRIC_VARS, +}; + +/* indices of the interpolated values of the above grid functions and their derivatives */ +enum InterpMetricVars { + I_GTXX = 0, + I_GTYY, + I_GTZZ, + I_GTXY, + I_GTXZ, + I_GTYZ, + I_PHI, + I_PHI_DX, + I_PHI_DY, + I_PHI_DZ, + I_ATXX, + I_ATYY, + I_ATZZ, + I_ATXY, + I_ATXZ, + I_ATYZ, + I_K, + I_K_DX, + I_K_DY, + I_K_DZ, + I_XTX, + I_XTY, + I_XTZ, + I_BETAX, + I_BETAY, + I_BETAZ, + I_ALPHA, + I_ALPHA_DX, + I_ALPHA_DY, + I_ALPHA_DZ, + I_ALPHA_DXX, + I_ALPHA_DYY, + I_ALPHA_DZZ, + I_ALPHA_DXY, + I_ALPHA_DXZ, + I_ALPHA_DYZ, + I_KDOT_XX, + I_KDOT_YY, + I_KDOT_ZZ, + I_KDOT_XY, + I_KDOT_XZ, + I_KDOT_YZ, + I_XTDOT_X, + I_XTDOT_Y, + I_XTDOT_Z, + I_PHIDOT, + I_PHIDOT_DX, + I_PHIDOT_DY, + I_PHIDOT_DZ, + NB_INTERP_VARS, +}; + +/* a set of basis functions */ +typedef struct BasisSet { + /* evaluate the idx-th basis function at the specified point*/ + double (*eval) (double coord, int idx); + /* evaluate the first derivative of the idx-th basis function at the specified point*/ + double (*eval_diff1)(double coord, int idx); + /* evaluate the second derivative of the idx-th basis function at the specified point*/ + double (*eval_diff2)(double coord, int idx); + /** + * Get the idx-th collocation point for the specified order. + * idx runs from 0 to order - 1 (inclusive) + */ + double (*colloc_point)(int order, int idx); +} BasisSet; + +extern const BasisSet qms_cheb_basis; +extern const BasisSet qms_cheb_even_basis; +extern const BasisSet qms_full_basis; +extern const BasisSet qms_tb_even_basis; +extern const BasisSet qms_sb_even_basis; +extern const BasisSet qms_tl_basis; +extern const BasisSet qms_cos_even_basis; + +extern double scale_factor; + +typedef struct MGContext MGContext; +typedef struct SORContext SORContext; + +/* precomputed values for a given refined grid */ +typedef struct CoordPatch { + CCTK_REAL origin[3]; + CCTK_INT delta[3]; + CCTK_INT size[3]; + + // basis values on the grid + double *basis_val_r; + double *basis_val_z; + + double *transform_z; + double *transform_matrix; + double *transform_matrix1; + double *transform_tmp; + double *one; + + int y_idx; + + MGContext *mg; + SORContext *sor; +} CoordPatch; + +/* state and scratch storage for the BiCGSTAB solver */ +typedef struct BiCGSTABContext { + double *p, *v, *y, *z, *t; + double *res, *res0; + double *k; + + cl_mem cl_p, cl_v, cl_y, cl_z, cl_t; + cl_mem cl_res, cl_res0; + cl_mem cl_k, cl_mat; + cl_mem cl_rho, cl_alpha, cl_beta, cl_omega, cl_omega1; + cl_mem cl_tmp, cl_tmp1; + +} BiCGSTABContext; + +typedef struct MaximalSlicingContext { + cGH *gh; + const BasisSet *basis; + const BasisSet *basis1; + + struct { + double time; + double *coeffs; + } solution_cache[8]; + int nb_solutions; + + double *coeffs_eval; + + BiCGSTABContext bicgstab; + int steps_since_inverse; + + uint64_t solve_count; + uint64_t solve_time; + + uint64_t lu_solves_count; + uint64_t lu_solves_time; + + uint64_t cg_solve_count; + uint64_t cg_iter_count; + uint64_t cg_time_total; + + uint64_t interp_geometry_count; + uint64_t interp_geometry_time; + + uint64_t calc_eq_coeffs_count; + uint64_t calc_eq_coeffs_time; + + uint64_t construct_matrix_count; + uint64_t construct_matrix_time; + + uint64_t grid_expand_count; + uint64_t grid_expand_time; + + // the grid of collocation points + double *colloc_grid[2]; + + // interpolation parameters + int coord_system; + int interp_operator; + int interp_params; + + CCTK_REAL *interp_coords[3]; + + int interp_vars_indices[NB_METRIC_VARS]; + + CCTK_REAL *interp_values[NB_INTERP_VARS]; + CCTK_REAL *interp_values_prefilter[NB_INTERP_VARS]; + CCTK_INT interp_value_codes[NB_INTERP_VARS]; + + CCTK_REAL *metric_u[6]; + + CCTK_REAL *kij_kij; + CCTK_REAL *trk; + + int nb_coeffs_x; + int nb_coeffs_z; + int nb_coeffs; + + int nb_colloc_points_x; + int nb_colloc_points_z; + int nb_colloc_points; + + int colloc_grid_order_x; + int colloc_grid_order_z; + + double *mat; + double *mat_f; + double *rhs; + double *coeffs; + int *ipiv; + double *basis_x_val; + double *basis_x_dval; + double *basis_x_d2val; + + double *basis_z_val; + double *basis_z_dval; + double *basis_z_d2val; + + double *basis_val_00; + double *basis_val_20; + double *basis_val_02; + double *basis_val_11; + double *basis_val_10; + double *basis_val_01; + + double *eq_coeff_00; + double *eq_coeff_20; + double *eq_coeff_02; + double *eq_coeff_11; + double *eq_coeff_10; + double *eq_coeff_01; + + double *coeff_scale; + + double *input_filter; + + CoordPatch *patches; + int nb_patches; + + // OpenCL / CLBLAS stuff + cl_context cl_ctx; + cl_command_queue cl_queue; + cl_device_id ocl_device; + + cl_mem ocl_coeffs; +} MaximalSlicingContext; + +int qms_maximal_solve(MaximalSlicingContext *ms); + +#include <sys/time.h> +static inline int64_t gettime(void) +{ + struct timeval tv; + gettimeofday(&tv, NULL); + return (int64_t)tv.tv_sec * 1000000 + tv.tv_usec; +} + |