/* * Quasimaximal slicing -- actual solver code * Copyright (C) 2016 Anton Khirnov * * This program 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. * * This program 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 this program. If not, see . */ #include "common.h" #include #include #include #include #include #if HAVE_OPENCL #include #include #endif #include "cctk.h" #include "cctk_Timers.h" #include "util_Table.h" #include "basis.h" #include "pssolve.h" #include "qms_solve.h" #define NB_COEFFS(qms) (qms->nb_coeffs[0] * qms->nb_coeffs[1]) #define NB_COLLOC_POINTS(qms) (qms->nb_colloc_points[0] * qms->nb_colloc_points[1]) /* 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_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, }; struct QMSSolverPriv { PSSolveContext *ps_ctx; cGH *gh; int colloc_grid_order[2]; double *eq_coeffs[PSSOLVE_DIFF_ORDER_NB]; double *rhs; double *coeff_scale; // 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_INT interp_value_codes[NB_INTERP_VARS]; #if HAVE_OPENCL // OpenCL / CLBLAS stuff cl_context ocl_ctx; cl_command_queue ocl_queue; #endif uint64_t solve_count; uint64_t solve_time; uint64_t interp_geometry_count; uint64_t interp_geometry_time; uint64_t calc_eq_coeffs_count; uint64_t calc_eq_coeffs_time; }; /* mapping between our indices and thorn names */ static const char *metric_vars[] = { #if QMS_CCZ4 [GTXX] = "ML_CCZ4::gt11", [GTYY] = "ML_CCZ4::gt22", [GTZZ] = "ML_CCZ4::gt33", [GTXY] = "ML_CCZ4::gt12", [GTXZ] = "ML_CCZ4::gt13", [GTYZ] = "ML_CCZ4::gt23", [ATXX] = "ML_CCZ4::At11", [ATYY] = "ML_CCZ4::At22", [ATZZ] = "ML_CCZ4::At33", [ATXY] = "ML_CCZ4::At12", [ATXZ] = "ML_CCZ4::At13", [ATYZ] = "ML_CCZ4::At23", [PHI] = "ML_CCZ4::phi", [K] = "ML_CCZ4::trK", [XTX] = "ML_CCZ4::Xt1", [XTY] = "ML_CCZ4::Xt2", [XTZ] = "ML_CCZ4::Xt3", [BETAX] = "ML_CCZ4::beta1", [BETAY] = "ML_CCZ4::beta2", [BETAZ] = "ML_CCZ4::beta3", [ALPHA] = "ML_CCZ4::alpha", [KDOT_XX] = "ML_CCZ4::Kdot11", [KDOT_YY] = "ML_CCZ4::Kdot22", [KDOT_ZZ] = "ML_CCZ4::Kdot33", [KDOT_XY] = "ML_CCZ4::Kdot12", [KDOT_XZ] = "ML_CCZ4::Kdot13", [KDOT_YZ] = "ML_CCZ4::Kdot23", [XTDOT_X] = "ML_CCZ4::Xtdot1", [XTDOT_Y] = "ML_CCZ4::Xtdot2", [XTDOT_Z] = "ML_CCZ4::Xtdot3", [PHIDOT] = "ML_CCZ4::phidot", #else [GTXX] = "ML_BSSN::gt11", [GTYY] = "ML_BSSN::gt22", [GTZZ] = "ML_BSSN::gt33", [GTXY] = "ML_BSSN::gt12", [GTXZ] = "ML_BSSN::gt13", [GTYZ] = "ML_BSSN::gt23", [ATXX] = "ML_BSSN::At11", [ATYY] = "ML_BSSN::At22", [ATZZ] = "ML_BSSN::At33", [ATXY] = "ML_BSSN::At12", [ATXZ] = "ML_BSSN::At13", [ATYZ] = "ML_BSSN::At23", [PHI] = "ML_BSSN::phi", [K] = "ML_BSSN::trK", [XTX] = "ML_BSSN::Xt1", [XTY] = "ML_BSSN::Xt2", [XTZ] = "ML_BSSN::Xt3", [BETAX] = "ML_BSSN::beta1", [BETAY] = "ML_BSSN::beta2", [BETAZ] = "ML_BSSN::beta3", [ALPHA] = "ML_BSSN::alpha", //[ALPHA] = "ADMBase::alp", [KDOT_XX] = "ML_BSSN::Kdot11", [KDOT_YY] = "ML_BSSN::Kdot22", [KDOT_ZZ] = "ML_BSSN::Kdot33", [KDOT_XY] = "ML_BSSN::Kdot12", [KDOT_XZ] = "ML_BSSN::Kdot13", [KDOT_YZ] = "ML_BSSN::Kdot23", [XTDOT_X] = "ML_BSSN::Xtdot1", [XTDOT_Y] = "ML_BSSN::Xtdot2", [XTDOT_Z] = "ML_BSSN::Xtdot3", [PHIDOT] = "ML_BSSN::phidot", #endif }; /* mapping between the cactus grid values and interpolated values */ static const CCTK_INT interp_operation_indices[] = { [I_GTXX] = GTXX, [I_GTYY] = GTYY, [I_GTZZ] = GTZZ, [I_GTXY] = GTXY, [I_GTXZ] = GTXZ, [I_GTYZ] = GTYZ, [I_PHI] = PHI, [I_PHI_DX] = PHI, [I_PHI_DY] = PHI, [I_PHI_DZ] = PHI, [I_ATXX] = ATXX, [I_ATYY] = ATYY, [I_ATZZ] = ATZZ, [I_ATXY] = ATXY, [I_ATXZ] = ATXZ, [I_ATYZ] = ATYZ, [I_K] = K, [I_XTX] = XTX, [I_XTY] = XTY, [I_XTZ] = XTZ, [I_BETAX] = BETAX, [I_BETAY] = BETAY, [I_BETAZ] = BETAZ, [I_ALPHA] = ALPHA, [I_ALPHA_DX] = ALPHA, [I_ALPHA_DY] = ALPHA, [I_ALPHA_DZ] = ALPHA, [I_ALPHA_DXX] = ALPHA, [I_ALPHA_DYY] = ALPHA, [I_ALPHA_DZZ] = ALPHA, [I_ALPHA_DXY] = ALPHA, [I_ALPHA_DXZ] = ALPHA, [I_ALPHA_DYZ] = ALPHA, [I_KDOT_XX] = KDOT_XX, [I_KDOT_YY] = KDOT_YY, [I_KDOT_ZZ] = KDOT_ZZ, [I_KDOT_XY] = KDOT_XY, [I_KDOT_XZ] = KDOT_XZ, [I_KDOT_YZ] = KDOT_YZ, [I_XTDOT_X] = XTDOT_X, [I_XTDOT_Y] = XTDOT_Y, [I_XTDOT_Z] = XTDOT_Z, [I_PHIDOT] = PHIDOT, [I_PHIDOT_DX] = PHIDOT, [I_PHIDOT_DY] = PHIDOT, [I_PHIDOT_DZ] = PHIDOT, }; /* the operation (plain value or x/y/z-derivative) to apply during interpolation */ static const CCTK_INT interp_operation_codes[] = { [I_GTXX] = 0, [I_GTYY] = 0, [I_GTZZ] = 0, [I_GTXY] = 0, [I_GTXZ] = 0, [I_GTYZ] = 0, [I_PHI] = 0, [I_PHI_DX] = 1, [I_PHI_DY] = 2, [I_PHI_DZ] = 3, [I_ATXX] = 0, [I_ATYY] = 0, [I_ATZZ] = 0, [I_ATXY] = 0, [I_ATXZ] = 0, [I_ATYZ] = 0, [I_K] = 0, [I_XTX] = 0, [I_XTY] = 0, [I_XTZ] = 0, [I_BETAX] = 0, [I_BETAY] = 0, [I_BETAZ] = 0, [I_ALPHA] = 0, [I_ALPHA_DX] = 1, [I_ALPHA_DY] = 2, [I_ALPHA_DZ] = 3, [I_ALPHA_DXX] = 11, [I_ALPHA_DYY] = 22, [I_ALPHA_DZZ] = 33, [I_ALPHA_DXY] = 12, [I_ALPHA_DXZ] = 13, [I_ALPHA_DYZ] = 23, [I_KDOT_XX] = 0, [I_KDOT_YY] = 0, [I_KDOT_ZZ] = 0, [I_KDOT_XY] = 0, [I_KDOT_XZ] = 0, [I_KDOT_YZ] = 0, [I_XTDOT_X] = 0, [I_XTDOT_Y] = 0, [I_XTDOT_Z] = 0, [I_PHIDOT] = 0, [I_PHIDOT_DX] = 1, [I_PHIDOT_DY] = 2, [I_PHIDOT_DZ] = 3, }; /* interpolate the cactus gridfunctions onto the pseudospectral grid */ static int interp_geometry(QMSSolver *ctx) { QMSSolverPriv *s = ctx->priv; int ret; ret = CCTK_InterpGridArrays(s->gh, 3, s->interp_operator, s->interp_params, s->coord_system, NB_COLLOC_POINTS(ctx), CCTK_VARIABLE_REAL, (const void * const *)s->interp_coords, ARRAY_ELEMS(s->interp_vars_indices), s->interp_vars_indices, ARRAY_ELEMS(s->interp_values), s->interp_value_codes, (void * const *)s->interp_values); if (ret < 0) CCTK_WARN(0, "Error interpolating"); return 0; } /* evaluate the equation coefficients at the collocation points */ static int calc_eq_coeffs(QMSSolver *ctx) { QMSSolverPriv *s = ctx->priv; //#pragma omp parallel for schedule(dynamic, ms->nb_colloc_points_x) for (int i = 0; i < NB_COLLOC_POINTS(ctx); i++) { const double x = s->interp_coords[0][i]; const double z = s->interp_coords[2][i]; const int zaxis = x <= EPS; double Am[3][3], K[3][3], Km[3][3], Ku[3][3], gtu[3][3]; double k2, kij_dij_alpha, k_kdot, k3; const double gtxx = s->interp_values[I_GTXX][i]; const double gtyy = s->interp_values[I_GTYY][i]; const double gtzz = s->interp_values[I_GTZZ][i]; const double gtxy = s->interp_values[I_GTXY][i]; const double gtxz = s->interp_values[I_GTXZ][i]; const double gtyz = s->interp_values[I_GTYZ][i]; const double gt[3][3] = {{ gtxx, gtxy, gtxz }, { gtxy, gtyy, gtyz }, { gtxz, gtyz, gtzz }}; const double Atxx = s->interp_values[I_ATXX][i]; const double Atyy = s->interp_values[I_ATYY][i]; const double Atzz = s->interp_values[I_ATZZ][i]; const double Atxy = s->interp_values[I_ATXY][i]; const double Atxz = s->interp_values[I_ATXZ][i]; const double Atyz = s->interp_values[I_ATYZ][i]; const double phi = s->interp_values[I_PHI][i]; const double phidot = s->interp_values[I_PHIDOT][i]; const double phidot_dx = s->interp_values[I_PHIDOT_DX][i]; const double phidot_dz = s->interp_values[I_PHIDOT_DZ][i]; const double At[3][3] = {{ Atxx, Atxy, Atxz }, { Atxy, Atyy, Atyz }, { Atxz, Atyz, Atzz }}; const double trK = s->interp_values[I_K][i]; const double kdot_xx = s->interp_values[I_KDOT_XX][i]; const double kdot_yy = s->interp_values[I_KDOT_YY][i]; const double kdot_zz = s->interp_values[I_KDOT_ZZ][i]; const double kdot_xy = s->interp_values[I_KDOT_XY][i]; const double kdot_xz = s->interp_values[I_KDOT_XZ][i]; const double kdot_yz = s->interp_values[I_KDOT_YZ][i]; const double kdot[3][3] = {{ kdot_xx, kdot_xy, kdot_xz }, { kdot_xy, kdot_yy, kdot_yz }, { kdot_xz, kdot_yz, kdot_zz }}; const double alpha = s->interp_values[I_ALPHA][i]; const double dx_alpha = s->interp_values[I_ALPHA_DX][i]; const double dz_alpha = s->interp_values[I_ALPHA_DZ][i]; const double dxx_alpha = s->interp_values[I_ALPHA_DXX][i]; const double dzz_alpha = s->interp_values[I_ALPHA_DZZ][i]; const double dxz_alpha = s->interp_values[I_ALPHA_DXZ][i]; const double dij_alpha[3][3] = {{ dxx_alpha, 0, dxz_alpha }, { 0, zaxis ? dxx_alpha : dx_alpha / x, 0 }, { dxz_alpha, 0, dzz_alpha }}; const double Xtx = s->interp_values[I_XTX][i]; const double Xtz = s->interp_values[I_XTZ][i]; const double Xtdot_x = s->interp_values[I_XTDOT_X][i]; const double Xtdot_z = s->interp_values[I_XTDOT_Z][i]; const double det = gtxx * gtyy * gtzz + 2 * gtxy * gtyz * gtxz - gtzz * SQR(gtxy) - SQR(gtxz) * gtyy - gtxx * SQR(gtyz); // \tilde{γ}^{ij} gtu[0][0] = (gtyy * gtzz - SQR(gtyz)) / det; gtu[1][1] = (gtxx * gtzz - SQR(gtxz)) / det; gtu[2][2] = (gtxx * gtyy - SQR(gtxy)) / det; gtu[0][1] = -(gtxy * gtzz - gtyz * gtxz) / det; gtu[0][2] = (gtxy * gtyz - gtyy * gtxz) / det; gtu[1][2] = -(gtxx * gtyz - gtxy * gtxz) / det; gtu[1][0] = gtu[0][1]; gtu[2][0] = gtu[0][2]; gtu[2][1] = gtu[1][2]; // K_{ij} for (int j = 0; j < 3; j++) for (int k = 0; k < 3; k++) K[j][k] = At[j][k] / SQR(phi) + gt[j][k] * trK; for (int j = 0; j < 3; j++) for (int k = 0; k < 3; k++) { double val = 0.0; for (int l = 0; l < 3; l++) val += SQR(phi) * gtu[j][l] * K[l][k]; Km[j][k] = val; } // K^{ij} for (int j = 0; j < 3; j++) for (int k = 0; k < 3; k++) { double val = 0.0; for (int l = 0; l < 3; l++) val += SQR(phi) * gtu[j][l] * Km[k][l]; Ku[j][k] = val; } // \tilde{A}_{i}^j for (int j = 0; j < 3; j++) for (int k = 0; k < 3; k++) { double val = 0.0; for (int l = 0; l < 3; l++) val += gtu[j][l] * At[l][k]; Am[j][k] = val; } kij_dij_alpha = 0.0; for (int j = 0; j < 3; j++) for (int k = 0; k < 3; k++) kij_dij_alpha += Ku[j][k] * dij_alpha[j][k]; k_kdot = 0.0; for (int j = 0; j < 3; j++) for (int k = 0; k < 3; k++) k_kdot += kdot[j][k] * Ku[j][k]; k3 = 0.0; for (int j = 0; j < 3; j++) for (int k = 0; k < 3; k++) { double val = 0.0; for (int l = 0; l < 3; l++) val += Km[k][l] * Ku[l][j]; k3 += val * K[j][k]; } // K_{ij} K^{ij} k2 = 0.0; for (int j = 0; j < 3; j++) for (int k = 0; k < 3; k++) k2 += Km[j][k] * Km[k][j]; { const double gtuxx = gtu[0][0]; const double gtuyy = gtu[1][1]; const double gtuzz = gtu[2][2]; const double gtuxz = gtu[0][2]; const double phi_dx = s->interp_values[I_PHI_DX][i]; const double phi_dz = s->interp_values[I_PHI_DZ][i]; const double Xtx = s->interp_values[I_XTX][i]; const double Xtz = s->interp_values[I_XTZ][i]; const double betax = s->interp_values[I_BETAX][i]; const double betaz = s->interp_values[I_BETAZ][i]; const double Xx = SQR(phi) * (Xtx + (phi_dx * gtuxx + phi_dz * gtuxz) / phi); const double Xz = SQR(phi) * (Xtz + (phi_dx * gtuxz + phi_dz * gtuzz) / phi); const double Xdot_x = 2 * phi * phidot * Xtx + SQR(phi) * Xtdot_x + phi * (phidot_dx * gtuxx + phidot_dz * gtuxz) - phidot * (phi_dx * gtuxx + phi_dz * gtuxz) + 2 * alpha * (phi_dx * Ku[0][0] + phi_dz * Ku[0][2]) / phi; const double Xdot_z = 2 * phi * phidot * Xtz + SQR(phi) * Xtdot_z + phi * (phidot_dz * gtuzz + phidot_dx * gtuxz) - phidot * (phi_dz * gtuzz + phi_dx * gtuxz) + 2 * alpha * (phi_dz * Ku[2][2] + phi_dx * Ku[0][2]) / phi; s->eq_coeffs[PSSOLVE_DIFF_ORDER_20][i] = SQR(phi) * (gtuxx + ((x <= EPS) ? gtuyy : 0.0)); s->eq_coeffs[PSSOLVE_DIFF_ORDER_02][i] = SQR(phi) * gtuzz; s->eq_coeffs[PSSOLVE_DIFF_ORDER_11][i] = SQR(phi) * gtuxz * 2; s->eq_coeffs[PSSOLVE_DIFF_ORDER_10][i] = -Xx + ((x > EPS) ? SQR(phi) * gtuyy / x : 0.0); s->eq_coeffs[PSSOLVE_DIFF_ORDER_01][i] = -Xz; s->eq_coeffs[PSSOLVE_DIFF_ORDER_00][i] = -k2; s->rhs[i] = -2 * alpha * kij_dij_alpha + Xdot_x * dx_alpha + Xdot_z * dz_alpha + 2 * (k_kdot + 2 * alpha * k3) * alpha; } } return 0; } int qms_solver_solve(QMSSolver *ctx) { QMSSolverPriv *s = ctx->priv; int ret; int64_t start, totaltime_start; totaltime_start = gettime(); /* interpolate the metric values and construct the quantities we'll need */ CCTK_TimerStart("QuasiMaximalSlicing_interp_geometry"); start = gettime(); ret = interp_geometry(ctx); s->interp_geometry_time += gettime() - start; s->interp_geometry_count++; CCTK_TimerStop("QuasiMaximalSlicing_interp_geometry"); if (ret < 0) return ret; CCTK_TimerStart("QuasiMaximalSlicing_calc_eq_coeffs"); start = gettime(); ret = calc_eq_coeffs(ctx); s->calc_eq_coeffs_time += gettime() - start; s->calc_eq_coeffs_count++; CCTK_TimerStop("QuasiMaximalSlicing_calc_eq_coeffs"); if (ret < 0) return ret; ret = qms_pssolve_solve(s->ps_ctx, (const double * const *)s->eq_coeffs, s->rhs, ctx->coeffs); if (ret < 0) return ret; for (int i = 0; i < NB_COEFFS(ctx); i++) ctx->coeffs[i] *= s->coeff_scale[i]; s->solve_count++; s->solve_time += gettime() - totaltime_start; return 0; } void qms_solver_print_stats(QMSSolver *ctx) { QMSSolverPriv *s = ctx->priv; fprintf(stderr, "%g%% interpolate geometry: %lu, " "total time %g s, avg time per call %g ms\n", (double)s->interp_geometry_time * 100 / s->solve_time, s->interp_geometry_count, (double)s->interp_geometry_time / 1e6, (double)s->interp_geometry_time / s->interp_geometry_count / 1e3); fprintf(stderr, "%g%% calc equation coefficients: %lu, " "total time %g s, avg time per call %g ms\n", (double)s->calc_eq_coeffs_time * 100 / s->solve_time, s->calc_eq_coeffs_count, (double)s->calc_eq_coeffs_time / 1e6, (double)s->calc_eq_coeffs_time / s->calc_eq_coeffs_count / 1e3); fprintf(stderr, "%g%% pseudospectral matrix construction: %lu, " "total time %g s, avg time per call %g ms\n", (double)s->ps_ctx->construct_matrix_time * 100 / s->solve_time, s->ps_ctx->construct_matrix_count, (double)s->ps_ctx->construct_matrix_time / 1e6, (double)s->ps_ctx->construct_matrix_time / s->ps_ctx->construct_matrix_count / 1e3); fprintf(stderr, "%g%% BiCGSTAB %lu solves, " "%lu iterations, total time %g s, " "avg iterations per solve %g, avg time per solve %g ms, " "avg time per iteration %g ms\n", (double)s->ps_ctx->cg_time_total * 100 / s->solve_time, s->ps_ctx->cg_solve_count, s->ps_ctx->cg_iter_count, (double)s->ps_ctx->cg_time_total / 1e6, (double)s->ps_ctx->cg_iter_count / s->ps_ctx->cg_solve_count, (double)s->ps_ctx->cg_time_total / s->ps_ctx->cg_solve_count / 1e3, (double)s->ps_ctx->cg_time_total / s->ps_ctx->cg_iter_count / 1e3); fprintf(stderr, "%g%% LU %lu solves, total time %g s, avg time per solve %g ms\n", (double)s->ps_ctx->lu_solves_time * 100 / s->solve_time, s->ps_ctx->lu_solves_count, (double)s->ps_ctx->lu_solves_time / 1e6, (double)s->ps_ctx->lu_solves_time / s->ps_ctx->lu_solves_count / 1e3); } static void init_opencl(QMSSolver *ctx) #if HAVE_OPENCL { QMSSolverPriv *s = ctx->priv; int err, count; cl_platform_id platform; cl_context_properties props[3]; cl_device_id ocl_device; err = clGetPlatformIDs(1, &platform, &count); if (err != CL_SUCCESS || count < 1) { fprintf(stderr, "Could not get an OpenCL platform ID\n"); return; } err = clGetDeviceIDs(platform, CL_DEVICE_TYPE_GPU, 1, &ocl_device, &count); if (err != CL_SUCCESS || count < 1) { fprintf(stderr, "Could not get an OpenCL device ID\n"); return; } props[0] = CL_CONTEXT_PLATFORM; props[1] = (cl_context_properties)platform; props[2] = 0; s->ocl_ctx = clCreateContext(props, 1, &ocl_device, NULL, NULL, &err); if (err != CL_SUCCESS || !s->ocl_ctx) { fprintf(stderr, "Could not create an OpenCL context\n"); return; } s->ocl_queue = clCreateCommandQueue(s->ocl_ctx, ocl_device, 0, &err); if (err != CL_SUCCESS || !s->ocl_queue) { fprintf(stderr, "Could not create an OpenCL command queue: %d\n", err); goto fail; } err = clblasSetup(); if (err != CL_SUCCESS) { fprintf(stderr, "Error setting up clBLAS\n"); goto fail; } return; fail: if (s->ocl_queue) clReleaseCommandQueue(s->ocl_queue); s->ocl_queue = 0; if (s->ocl_ctx) clReleaseContext(s->ocl_ctx); s->ocl_ctx = 0; } #else { } #endif int qms_solver_init(QMSSolver **pctx, cGH *cctkGH, int basis_order_r, int basis_order_z, double sf, double filter_power, double input_filter_power) { QMSSolver *ctx; QMSSolverPriv *s; int ret; ctx = calloc(1, sizeof(*ctx)); if (!ctx) return -ENOMEM; ctx->priv = calloc(1, sizeof(*ctx->priv)); if (!ctx->priv) goto fail; s = ctx->priv; s->gh = cctkGH; ctx->basis[0] = &qms_sb_even_basis; #if QMS_POLAR ctx->basis[1] = &qms_cos_even_basis; #else ctx->basis[1] = &qms_sb_even_basis; #endif ctx->nb_coeffs[0] = basis_order_r; ctx->nb_coeffs[1] = basis_order_z; ctx->nb_colloc_points[0] = basis_order_r; ctx->nb_colloc_points[1] = basis_order_z; if (NB_COLLOC_POINTS(ctx) != NB_COEFFS(ctx)) CCTK_WARN(0, "Non-square collocation matrix"); s->colloc_grid_order[0] = ctx->nb_colloc_points[0]; s->colloc_grid_order[1] = ctx->nb_colloc_points[1]; ret = posix_memalign((void**)&ctx->coeffs, 32, sizeof(*ctx->coeffs) * NB_COEFFS(ctx)); ret |= posix_memalign((void**)&s->rhs, 32, sizeof(*s->rhs) * NB_COLLOC_POINTS(ctx)); if (ret) goto fail; //FIXME scale_factor = 1.0; scale_factor = (64.0 / ctx->basis[0]->colloc_point(s->colloc_grid_order[0], ctx->nb_colloc_points[0] - 1)); fprintf(stderr, "scale factor %16.16g\n", scale_factor); init_opencl(ctx); ret = qms_pssolve_context_alloc(&s->ps_ctx); if (ret < 0) CCTK_WARN(0, "Error allocating the pseudospectral solver"); s->ps_ctx->basis[0] = ctx->basis[0]; s->ps_ctx->basis[1] = ctx->basis[1]; s->ps_ctx->solve_order[0] = basis_order_r; s->ps_ctx->solve_order[1] = basis_order_z; #if HAVE_OPENCL s->ps_ctx->ocl_ctx = s->ocl_ctx; s->ps_ctx->ocl_queue = s->ocl_queue; #endif ret = qms_pssolve_context_init(s->ps_ctx); if (ret < 0) CCTK_WARN(0, "Error initializing the pseudospectral solver"); for (int i = 0; i < MAX(s->ps_ctx->solve_order[0], s->ps_ctx->solve_order[1]); i++) { fprintf(stderr, "%d ", i); if (i < s->ps_ctx->solve_order[0]) fprintf(stderr, "%g\t", s->ps_ctx->colloc_grid[0][i]); else fprintf(stderr, "\t\t"); if (i < s->ps_ctx->solve_order[1]) fprintf(stderr, "%g\t", s->ps_ctx->colloc_grid[1][i]); fprintf(stderr, "\n"); } for (int i = 0; i < ARRAY_ELEMS(s->eq_coeffs); i++) { ret = posix_memalign((void**)&s->eq_coeffs[i], 32, NB_COLLOC_POINTS(ctx) * sizeof(*s->eq_coeffs[i])); if (ret) goto fail; } for (int i = 0; i < ARRAY_ELEMS(s->interp_coords); i++) { ret |= posix_memalign((void**)&s->interp_coords[i], 32, NB_COLLOC_POINTS(ctx) * sizeof(*s->interp_coords[i])); } if (ret) goto fail; for (int i = 0; i < ctx->nb_colloc_points[1]; i++) { for (int j = 0; j < ctx->nb_colloc_points[0]; j++) { #if QMS_POLAR double phi = s->ps_ctx->colloc_grid[1][i]; double r = s->ps_ctx->colloc_grid[0][j]; double x = r * cos(phi); double z = r * sin(phi); #else double x = s->ps_ctx->colloc_grid[0][j]; double z = s->ps_ctx->colloc_grid[1][i]; #endif s->interp_coords[0][i * ctx->nb_colloc_points[0] + j] = x; s->interp_coords[1][i * ctx->nb_colloc_points[0] + j] = 0; s->interp_coords[2][i * ctx->nb_colloc_points[0] + j] = z; } } ret = posix_memalign((void**)&s->coeff_scale, 32, NB_COEFFS(ctx) * sizeof(*s->coeff_scale)); if (ret) goto fail; for (int j = 0; j < ctx->nb_coeffs[1]; j++) for (int i = 0; i < ctx->nb_coeffs[0]; i++) { s->coeff_scale[j * ctx->nb_coeffs[0] + i] = exp(-36.0 * pow((double)i / ctx->nb_coeffs[0], filter_power)) * exp(-36.0 * pow((double)j / ctx->nb_coeffs[1], filter_power)); } for (int i = 0; i < ARRAY_ELEMS(s->interp_values); i++) { ret = posix_memalign((void**)&s->interp_values[i], 32, NB_COLLOC_POINTS(ctx) * sizeof(*s->interp_values[i])); if (ret) goto fail; s->interp_value_codes[i] = CCTK_VARIABLE_REAL; } for (int i = 0; i < ARRAY_ELEMS(metric_vars); i++) { s->interp_vars_indices[i] = CCTK_VarIndex(metric_vars[i]); if (s->interp_vars_indices[i] < 0) CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, "Error getting the index of variable: %s\n", metric_vars[i]); } s->coord_system = CCTK_CoordSystemHandle("cart3d"); if (s->coord_system < 0) CCTK_WARN(0, "Error getting the coordinate system"); s->interp_operator = CCTK_InterpHandle("Lagrange polynomial interpolation (tensor product)"); if (s->interp_operator < 0) CCTK_WARN(0, "Error getting the interpolation operator"); s->interp_params = Util_TableCreateFromString("order=4 want_global_mode=1"); if (s->interp_params < 0) CCTK_WARN(0, "Error creating interpolation parameters table"); ret = Util_TableSetIntArray(s->interp_params, NB_INTERP_VARS, interp_operation_codes, "operation_codes"); if (ret < 0) CCTK_WARN(0, "Error setting operation codes"); ret = Util_TableSetIntArray(s->interp_params, NB_INTERP_VARS, interp_operation_indices, "operand_indices"); if (ret < 0) CCTK_WARN(0, "Error setting operand indices"); CCTK_TimerCreate("QuasiMaximalSlicing_Solve"); CCTK_TimerCreate("QuasiMaximalSlicing_Expand"); CCTK_TimerCreate("QuasiMaximalSlicing_interp_geometry"); CCTK_TimerCreate("QuasiMaximalSlicing_calc_eq_coeffs"); CCTK_TimerCreate("QuasiMaximalSlicing_construct_matrix"); CCTK_TimerCreate("QuasiMaximalSlicing_solve_LU"); CCTK_TimerCreate("QuasiMaximalSlicing_solve_BiCGSTAB"); *pctx = ctx; return 0; fail: qms_solver_free(&ctx); return -ENOMEM; } void qms_solver_free(QMSSolver **pctx) { QMSSolver *ctx = *pctx; if (!ctx) return; if (ctx->priv) { for (int i = 0; i < ARRAY_ELEMS(ctx->priv->interp_coords); i++) free(ctx->priv->interp_coords[i]); for (int i = 0; i < ARRAY_ELEMS(ctx->priv->interp_values); i++) free(ctx->priv->interp_values[i]); for (int i = 0; i < ARRAY_ELEMS(ctx->priv->eq_coeffs); i++) free(ctx->priv->eq_coeffs[i]); free(ctx->priv->rhs); free(ctx->priv->coeff_scale); qms_pssolve_context_free(&ctx->priv->ps_ctx); #if HAVE_OPENCL if (ctx->priv->ocl_queue) clReleaseCommandQueue(ctx->priv->ocl_queue); if (ctx->priv->ocl_ctx) clReleaseContext(ctx->priv->ocl_ctx); #endif } free(ctx->priv); free(ctx->coeffs); free(ctx); *pctx = NULL; }