path: root/src/maximal_slicing_axi_mg.c

#include <ctype.h>
#include <errno.h>
#include <float.h>
#include <inttypes.h>
#include <limits.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include <mg2d.h>

#include <cblas.h>

#include "cctk.h"
#include "cctk_Arguments.h"
#include "cctk_Parameters.h"
#include "cctk_Timers.h"
#include "util_Table.h"

#define SQR(x) ((x) * (x))
#define ARRAY_ELEMS(x) (sizeof(x) / sizeof(*x))

/* precomputed values for a given refined grid */
typedef struct CoordPatch {
    double origin[3];
    int    delta[3];
    int    size[3];
    int    level;

    MG2DContext *solver;

    ptrdiff_t y_idx;

    /* number of x/z grid points by which the elliptic solver domain is offset
     * from the cactus grid */
    ptrdiff_t offset_left[2];
    ptrdiff_t offset_right[2];

    /* boundary interpolator parameters */
    double *interp_coords[3];

    double *boundary_interp_vals;
    size_t nb_boundary_interp_vals;

    int interp_values_codes[1];
    int interp_var_indices[1];
    int interp_operation_codes[1];
    int interp_operation_indices[1];

    int coord_system;
    int interp_operator;
    int interp_params;

    /**
     * pointers into boundary_interp_vals
     * the first row are the boundary points for the solver
     * [0] is the upper-z patch
     *     rows run from x=0 to x[lsh[0] - 1]
     * [1] is the upper-x patch
     *     rows run from z=0 to z=<upper physical boundary>
     */
    double   *boundary_val[2];
    ptrdiff_t boundary_val_stride[2];

} CoordPatch;

typedef struct MSMGContext {
    cGH *gh;

    CoordPatch *patches;
    int nb_patches;
} MSMGContext;

static int ctz(int a)
{
    int ret = 0;

    if (!a)
        return INT_MAX;

    while (!(a & 1)) {
        a >>= 1;
        ret++;
    }

    return ret;
}

static void coord_patch_free(CoordPatch *cp)
{
    mg2d_solver_free(&cp->solver);
}

static CoordPatch *get_coord_patch(MSMGContext *ms)
{
    cGH *gh = ms->gh;

    const double *a_x = CCTK_VarDataPtr(gh, 0, "grid::x");
    const double *a_y = CCTK_VarDataPtr(gh, 0, "grid::y");
    const double *a_z = CCTK_VarDataPtr(gh, 0, "grid::z");

    CoordPatch *cp;
    size_t domain_size[2];
    int integrator_substeps = MoLNumIntegratorSubsteps();
    int i, ret;

    for (int i = 0; i < ms->nb_patches; i++) {
        cp = &ms->patches[i];

        if (cp->origin[0] == ms->gh->cctk_origin_space[0] &&
            cp->origin[1] == ms->gh->cctk_origin_space[1] &&
            cp->origin[2] == ms->gh->cctk_origin_space[2] &&
            cp->size[0]   == ms->gh->cctk_lsh[0]          &&
            cp->size[1]   == ms->gh->cctk_lsh[1]          &&
            cp->size[2]   == ms->gh->cctk_lsh[2]          &&
            cp->delta[0]  == ms->gh->cctk_levfac[0]  &&
            cp->delta[1]  == ms->gh->cctk_levfac[1]  &&
            cp->delta[2]  == ms->gh->cctk_levfac[2])
            return cp;
    }

    /* create a new patch */
    ms->patches = realloc(ms->patches, sizeof(*ms->patches) * (ms->nb_patches + 1));
    cp = &ms->patches[ms->nb_patches];

    memset(cp, 0, sizeof(*cp));

    memcpy(cp->origin, ms->gh->cctk_origin_space, sizeof(cp->origin));
    memcpy(cp->size,   ms->gh->cctk_lsh,          sizeof(cp->size));
    memcpy(cp->delta,  ms->gh->cctk_levfac,  sizeof(cp->delta));
    cp->level = ctz(ms->gh->cctk_levfac[0]);

    for (i = 0; i < cp->size[1]; i++)
        if (fabs(a_y[CCTK_GFINDEX3D(gh, 0, i, 0)]) < 1e-8) {
            cp->y_idx = i;
            break;
        }
    if (i == cp->size[1])
        CCTK_WARN(0, "The grid does not include y==0");

    for (int i = 0; i < 2; i++) {
        ptrdiff_t offset_left = ms->gh->cctk_nghostzones[2 * i];
        ptrdiff_t offset_right = offset_left;

        /* account for grid tapering on refined levels */
        if (cp->level > 0)
            offset_right *= integrator_substeps * 2;

        domain_size[i] = ms->gh->cctk_lsh[2 * i] - offset_left - offset_right;
        /* the outer boundary layer overlaps with the first ghost zone */
        if (cp->level > 0)
            domain_size[i]++;

        cp->offset_left[i]  = offset_left;
        cp->offset_right[i] = offset_right;
    }

    if (domain_size[0] != domain_size[1])
        CCTK_WARN(0, "The grid is non-square, only square grids are supported");

    cp->solver = mg2d_solver_alloc(domain_size[0]);
    if (!cp->solver)
        CCTK_WARN(0, "Error allocating the solver");

    cp->solver->step[0] = a_x[1] - a_x[0];
    cp->solver->step[1] = cp->solver->step[0];

    cp->solver->fd_stencil = 1;

    cp->solver->boundaries[MG2D_BOUNDARY_0L]->type = MG2D_BC_TYPE_FIXDIFF;
    cp->solver->boundaries[MG2D_BOUNDARY_1L]->type = MG2D_BC_TYPE_FIXDIFF;
    cp->solver->boundaries[MG2D_BOUNDARY_0U]->type = MG2D_BC_TYPE_FIXVAL;
    cp->solver->boundaries[MG2D_BOUNDARY_1U]->type = MG2D_BC_TYPE_FIXVAL;

    cp->solver->maxiter   = 32;
    cp->solver->tol       = 1e-10;
    cp->solver->nb_cycles = 2;

    /* initialize boundary values to zero,
     * non-zero values on the outer boundaries of refined levels are filled in elsewhere */
    for (int i = 0; i < 4; i++)
        memset(cp->solver->boundaries[i]->val, 0, sizeof(double) * cp->solver->boundaries[i]->val_len);

    /* initialize the interpolation state */
    if (cp->level > 0) {
        cp->boundary_val_stride[0] = gh->cctk_lsh[0] - cp->offset_left[0];
        cp->boundary_val_stride[1] = cp->solver->domain_size;

        cp->nb_boundary_interp_vals = cp->boundary_val_stride[0] * cp->offset_right[1] +
                                      cp->boundary_val_stride[1] * cp->offset_right[0];

        ret = posix_memalign((void**)&cp->boundary_interp_vals, 32,
                             cp->nb_boundary_interp_vals * sizeof(*cp->boundary_interp_vals));
        if (ret != 0)
            CCTK_WARN(0, "Error allocating arrays for boundary values");
        cp->boundary_val[0] = cp->boundary_interp_vals;
        cp->boundary_val[1] = cp->boundary_interp_vals + cp->boundary_val_stride[0] * cp->offset_right[1];

        for (int i = 0; i < 3; i++) {
            ret = posix_memalign((void**)&cp->interp_coords[i], 32,
                                 cp->nb_boundary_interp_vals * sizeof(*cp->interp_coords[i]));
            if (ret != 0)
                CCTK_WARN(0, "Error allocating interp coords");
        }

        for (ptrdiff_t idx_row = 0; idx_row < cp->offset_right[1]; idx_row++)
            for (ptrdiff_t i = 0; i < cp->boundary_val_stride[0]; i++) {
                ptrdiff_t idx_dst = idx_row * cp->boundary_val_stride[0] + i;
                ptrdiff_t idx_src_x = CCTK_GFINDEX3D(gh, cp->offset_left[0] + i, 0, 0);
                ptrdiff_t idx_src_z = CCTK_GFINDEX3D(gh, 0, 0, gh->cctk_lsh[2] - cp->offset_right[1] + idx_row);

                cp->interp_coords[0][idx_dst] = a_x[idx_src_x];
                cp->interp_coords[1][idx_dst] = 0.0;
                cp->interp_coords[2][idx_dst] = a_z[idx_src_z];
            }

        for (ptrdiff_t idx_row = 0; idx_row < cp->offset_right[0]; idx_row++)
            for (ptrdiff_t i = 0; i < cp->boundary_val_stride[1]; i++) {
                ptrdiff_t idx_dst = (cp->boundary_val[1] - cp->boundary_val[0]) +
                                    idx_row * cp->boundary_val_stride[1] + i;
                ptrdiff_t idx_src_x = CCTK_GFINDEX3D(gh, gh->cctk_lsh[0] - cp->offset_right[0] + idx_row,
                                                     0, 0);
                ptrdiff_t idx_src_z = CCTK_GFINDEX3D(gh, 0, 0, cp->offset_left[1] + i);

                cp->interp_coords[0][idx_dst] = a_x[idx_src_x];
                cp->interp_coords[1][idx_dst] = 0.0;
                cp->interp_coords[2][idx_dst] = a_z[idx_src_z];
            }

#if 0
        for (int i = 0; i < gh->cctk_lsh[0]; i++) {
            cp->interp_coords[0][i] = a_x[CCTK_GFINDEX3D(gh, i, 0, 0)];
            cp->interp_coords[1][i] = 0.0;
            cp->interp_coords[2][i] = a_z[CCTK_GFINDEX3D(gh, 0, 0, gh->cctk_lsh[2] - gh->cctk_nghostzones[2] - 1)];
        }
        for (int i = 0; i < gh->cctk_lsh[2]; i++) {
            int idx = gh->cctk_lsh[0] + i;
            cp->interp_coords[0][idx] = a_x[CCTK_GFINDEX3D(gh, gh->cctk_lsh[0] - gh->cctk_nghostzones[0] - 1, 0, 0)];
            cp->interp_coords[1][idx] = 0.0;
            cp->interp_coords[2][idx] = a_z[CCTK_GFINDEX3D(gh, 0, 0, i)];
        }
#endif

        cp->interp_values_codes[0] = CCTK_VARIABLE_REAL;
        cp->interp_var_indices[0]  = CCTK_VarIndex("ML_BSSN::alpha");
        if (cp->interp_var_indices[0] < 0)
            CCTK_WARN(0, "Error getting the index of lapse variable");
        cp->interp_operation_codes[0] = 0;
        cp->interp_operation_indices[0] = 0;

        cp->coord_system = CCTK_CoordSystemHandle("cart3d");
        if (cp->coord_system < 0)
            CCTK_WARN(0, "Error getting the coordinate system");

        cp->interp_operator = CCTK_InterpHandle("Lagrange polynomial interpolation (tensor product)");
        if (cp->interp_operator < 0)
            CCTK_WARN(0, "Error getting the interpolation operator");

        cp->interp_params = Util_TableCreateFromString("order=4 want_global_mode=1");
        if (cp->interp_params < 0)
            CCTK_WARN(0, "Error creating interpolation parameters table");

        ret = Util_TableSetInt(cp->interp_params, cp->level, "max_reflevel");
        if (ret < 0)
            CCTK_WARN(0, "Error setting maximum reflevel for interpolation");

        ret = Util_TableSetIntArray(cp->interp_params, ARRAY_ELEMS(cp->interp_operation_codes),
                                    cp->interp_operation_codes, "operation_codes");
        if (ret < 0)
            CCTK_WARN(0, "Error setting operation codes");

        ret = Util_TableSetIntArray(cp->interp_params, ARRAY_ELEMS(cp->interp_operation_indices),
                                    cp->interp_operation_indices, "operand_indices");
        if (ret < 0)
            CCTK_WARN(0, "Error setting operand indices");
    }

    ms->nb_patches++;
    return cp;
}

static int context_init(cGH *gh, MSMGContext **ctx)
{
    MSMGContext *ms;
    int ret;

    ms = calloc(1, sizeof(*ms));
    if (!ms)
        return -ENOMEM;

    ms->gh = gh;

    *ctx = ms;

    return 0;
}

static void context_free(MSMGContext **pms)
{
    MSMGContext *ms = *pms;

    if (!ms)
        return;

    for (int i = 0; i < ms->nb_patches; i++)
        coord_patch_free(&ms->patches[i]);
    free(ms->patches);

    free(ms);
    *pms = NULL;
}

static void fill_eq_coeffs(cGH *gh, CoordPatch *cp)
{
    int ret;

    const ptrdiff_t stride_z = CCTK_GFINDEX3D(gh, 0, 0, 1);

    double *a_x = CCTK_VarDataPtr(gh, 0, "grid::x");
    double *a_z = CCTK_VarDataPtr(gh, 0, "grid::z");

    const double dx = a_x[1] - a_x[0];
    const double dz = a_z[stride_z] - a_z[0];

    double *a_gtxx = CCTK_VarDataPtr(gh, 0, "ML_BSSN::gt11");
    double *a_gtyy = CCTK_VarDataPtr(gh, 0, "ML_BSSN::gt22");
    double *a_gtzz = CCTK_VarDataPtr(gh, 0, "ML_BSSN::gt33");
    double *a_gtxy = CCTK_VarDataPtr(gh, 0, "ML_BSSN::gt12");
    double *a_gtxz = CCTK_VarDataPtr(gh, 0, "ML_BSSN::gt13");
    double *a_gtyz = CCTK_VarDataPtr(gh, 0, "ML_BSSN::gt23");
    double *a_Atxx = CCTK_VarDataPtr(gh, 0, "ML_BSSN::At11");
    double *a_Atyy = CCTK_VarDataPtr(gh, 0, "ML_BSSN::At22");
    double *a_Atzz = CCTK_VarDataPtr(gh, 0, "ML_BSSN::At33");
    double *a_Atxy = CCTK_VarDataPtr(gh, 0, "ML_BSSN::At12");
    double *a_Atxz = CCTK_VarDataPtr(gh, 0, "ML_BSSN::At13");
    double *a_Atyz = CCTK_VarDataPtr(gh, 0, "ML_BSSN::At23");
    double  *a_phi = CCTK_VarDataPtr(gh, 0, "ML_BSSN::phi");
    double  *a_Xtx = CCTK_VarDataPtr(gh, 0, "ML_BSSN::Xt1");
    double  *a_Xtz = CCTK_VarDataPtr(gh, 0, "ML_BSSN::Xt3");

    for (int idx_z = 0; idx_z < cp->solver->domain_size; idx_z++)
        for (int idx_x = 0; idx_x < cp->solver->domain_size; idx_x++) {
            const int idx_src = CCTK_GFINDEX3D(gh, idx_x + gh->cctk_nghostzones[0], cp->y_idx, idx_z + gh->cctk_nghostzones[2]);
            const int idx_dc  = idx_z * cp->solver->diff_coeffs_stride + idx_x;
            const int idx_rhs = idx_z * cp->solver->rhs_stride + idx_x;

            const double x = a_x[idx_src];

            const double gtxx = a_gtxx[idx_src];
            const double gtyy = a_gtyy[idx_src];
            const double gtzz = a_gtzz[idx_src];
            const double gtxy = a_gtxy[idx_src];
            const double gtxz = a_gtxz[idx_src];
            const double gtyz = a_gtyz[idx_src];
            const double Atxx = a_Atxx[idx_src];
            const double Atyy = a_Atyy[idx_src];
            const double Atzz = a_Atzz[idx_src];
            const double Atxy = a_Atxy[idx_src];
            const double Atxz = a_Atxz[idx_src];
            const double Atyz = a_Atyz[idx_src];
            const double  phi =  a_phi[idx_src];
            const double  Xtx =  a_Xtx[idx_src];
            const double  Xtz =  a_Xtz[idx_src];

            const double phi_dx = (a_phi[idx_src + 1] - a_phi[idx_src - 1]) / (2.0 * dx);
            const double phi_dz = (a_phi[idx_src + stride_z] - a_phi[idx_src - stride_z]) / (2.0 * dz);

            const double det = gtxx * gtyy * gtzz + 2 * gtxy * gtyz * gtxz - gtzz * SQR(gtxy) - SQR(gtxz) * gtyy - gtxx * SQR(gtyz);

            const double At[3][3] = {
                { Atxx, Atxy, Atxz },
                { Atxy, Atyy, Atyz },
                { Atxz, Atyz, Atzz }};


            double Xx, Xz, k2;

            double Am[3][3], gtu[3][3];

			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];

			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;
				}

			// K_{ij} K^{ij}
			k2 = 0.0;
			for (int j = 0; j < 3; j++)
				for (int k = 0; k < 3; k++)
					k2 += Am[j][k] * Am[k][j];

            Xx = SQR(phi) * (Xtx + (phi_dx * gtu[0][0] + phi_dz * gtu[0][2]) / phi);
            Xz = SQR(phi) * (Xtz + (phi_dx * gtu[0][2] + phi_dz * gtu[2][2]) / phi);

            cp->solver->diff_coeffs[MG2D_DIFF_COEFF_20][idx_dc] = SQR(phi) * (gtu[0][0] + ((idx_x == 0) ? gtu[1][1] : 0.0));
            cp->solver->diff_coeffs[MG2D_DIFF_COEFF_02][idx_dc] = SQR(phi) * gtu[2][2];
            cp->solver->diff_coeffs[MG2D_DIFF_COEFF_11][idx_dc] = SQR(phi) * gtu[0][2] * 2;
            cp->solver->diff_coeffs[MG2D_DIFF_COEFF_10][idx_dc] = -Xx + (idx_x ? SQR(phi) * gtu[1][1] / x : 0.0);
            cp->solver->diff_coeffs[MG2D_DIFF_COEFF_01][idx_dc] = -Xz;
            cp->solver->diff_coeffs[MG2D_DIFF_COEFF_00][idx_dc] = -k2;
            cp->solver->rhs[idx_rhs]                            = k2;
        }

    if (cp->level > 0) {
        ret = CCTK_InterpGridArrays(gh, 3, cp->interp_operator, cp->interp_params,
                                    cp->coord_system, cp->nb_boundary_interp_vals, CCTK_VARIABLE_REAL,
                                    (const void * const *)cp->interp_coords, ARRAY_ELEMS(cp->interp_var_indices),
                                    cp->interp_var_indices, ARRAY_ELEMS(cp->interp_values_codes), cp->interp_values_codes,
                                    (void * const *)&cp->boundary_interp_vals);
        if (ret < 0)
            CCTK_WARN(0, "Error interpolating");

        for (size_t i = 0; i < cp->solver->boundaries[MG2D_BOUNDARY_0U]->val_len; i++)
            cp->solver->boundaries[MG2D_BOUNDARY_0U]->val[i] = cp->boundary_val[0][i] - 1.0;
        for (size_t i = 0; i < cp->solver->boundaries[MG2D_BOUNDARY_1U]->val_len; i++)
            cp->solver->boundaries[MG2D_BOUNDARY_1U]->val[i] = cp->boundary_val[1][i] - 1.0;
    }
}

void maximal_slicing_axi_mg(CCTK_ARGUMENTS)
{
    static MSMGContext *ms;

    CoordPatch *cp;

    DECLARE_CCTK_ARGUMENTS;
    DECLARE_CCTK_PARAMETERS;

    double *coeffs = NULL;
    int i, ret;

    /* on the first run, init the solver */
    if (!ms)
        context_init(cctkGH, &ms);

    fprintf(stderr, "ms mg solve: level %d time %g\n", ctz(ms->gh->cctk_levfac[0]), ms->gh->cctk_time);

    cp = get_coord_patch(ms);

    fill_eq_coeffs(ms->gh, cp);

    CCTK_TimerStart("MaximalSlicingAxiMG_Solve");
    ret = mg2d_solve(cp->solver);
    if (ret < 0)
        CCTK_WARN(0, "Error solving the maximal slicing equation");
    CCTK_TimerStop("MaximalSlicingAxiMG_Solve");


    CCTK_TimerStart("MaximalSlicingAxiMG_Copy");
//#pragma omp parallel for
    for (int j = 0; j < cp->solver->domain_size; j++)
        for (int i = 0; i < cp->solver->domain_size; i++) {
            const int idx_dst = CCTK_GFINDEX3D(ms->gh, i + cp->offset_left[0], cp->y_idx, j + cp->offset_left[1]);
            const int idx_src = j * cp->solver->u_stride + i;
            alpha[idx_dst] = 1.0 + cp->solver->u[idx_src];
        }
    if (cp->level > 0) {
        /* fill in the interpolated outer buffer points */
        for (int idx_z = 0; idx_z < cp->offset_right[1]; idx_z++)
            for (int idx_x = 0; idx_x < cp->boundary_val_stride[0]; idx_x++) {
                const ptrdiff_t idx_src = idx_z * cp->boundary_val_stride[0] + idx_x;
                const ptrdiff_t idx_dst = CCTK_GFINDEX3D(ms->gh, cp->offset_left[0] + idx_x, cp->y_idx, ms->gh->cctk_lsh[2] - cp->offset_right[1] + idx_z);
                alpha[idx_dst] = cp->boundary_val[0][idx_src];
            }
        for (int idx_z = 0; idx_z < cp->boundary_val_stride[1]; idx_z++)
            for (int idx_x = 0; idx_x < cp->offset_right[0]; idx_x++) {
                const ptrdiff_t idx_src = idx_x * cp->boundary_val_stride[1] + idx_z;
                const ptrdiff_t idx_dst = CCTK_GFINDEX3D(ms->gh, ms->gh->cctk_lsh[0] - cp->offset_right[0] + idx_x, cp->y_idx,
                                                         cp->offset_left[1] + idx_z);
                alpha[idx_dst] = cp->boundary_val[1][idx_src];
            }
    }

    /* fill in the axial symmetry ghostpoints by mirroring */
    for (int idx_z = cp->offset_left[1]; idx_z < ms->gh->cctk_lsh[2]; idx_z++)
        for (int idx_x = 0; idx_x < cp->offset_left[0]; idx_x++) {
            const ptrdiff_t idx_dst = CCTK_GFINDEX3D(ms->gh, idx_x, cp->y_idx, idx_z);
            const ptrdiff_t idx_src = CCTK_GFINDEX3D(ms->gh, 2 * cp->offset_left[0] - idx_x, cp->y_idx, idx_z);
            alpha[idx_dst] = alpha[idx_src];
        }
    for (int idx_z = 0; idx_z < cp->offset_left[1]; idx_z++)
        for (int idx_x = 0; idx_x < ms->gh->cctk_lsh[0]; idx_x++) {
            const ptrdiff_t idx_dst = CCTK_GFINDEX3D(ms->gh, idx_x, cp->y_idx, idx_z);
            const ptrdiff_t idx_src = CCTK_GFINDEX3D(ms->gh, idx_x, cp->y_idx, 2 * cp->offset_left[1] - idx_z);
            alpha[idx_dst] = alpha[idx_src];
        }
    CCTK_TimerStop("MaximalSlicingAxiMG_Copy");
}