path: root/src/qms_solve.c

/*
 * Quasimaximal slicing -- actual solver code
 * Copyright (C) 2016 Anton Khirnov <anton@khirnov.net>
 *
 * 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 <http://www.gnu.org/licenses/>.
 */

#include "common.h"

#include <errno.h>
#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>

#if HAVE_OPENCL
#include <cl.h>
#include <clBLAS.h>
#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;
}