mhd_stat.c

Go to the documentation of this file.

#include <stdlib.h>
#include "../ascot5.h"
#include "../print.h"
#include "../error.h"
#include "../boozer.h"
#include "../spline/interp.h"
#include "../B_field.h"
#include "../math.h"
#include "../mhd.h"
#include "mhd_stat.h"

int mhd_stat_init(mhd_stat_data* data, int nmode, int nrho,
                  real rhomin, real rhomax, int* moden, int* modem,
                  real* amplitude_nm, real* omega_nm, real* phase_nm,
                  real* alpha, real* phi) {
 
    int err = 0;
    data->n_modes = nmode;
    data->rho_min = rhomin;
    data->rho_max = rhomax;
    data->nmode = (int*) malloc(nmode * sizeof(int));
    data->mmode = (int*) malloc(nmode * sizeof(int));
    data->omega_nm = (real*) malloc(nmode * sizeof(real));
    data->phase_nm = (real*) malloc(nmode * sizeof(real));
    data->amplitude_nm = (real*) malloc(nmode * sizeof(real));
    data->phi_nm = (interp1D_data*) malloc(nmode * sizeof(interp1D_data));
    data->alpha_nm = (interp1D_data*) malloc(nmode * sizeof(interp1D_data));
    for(int i = 0; i < nmode; i++) {
        data->nmode[i] = moden[i];
        data->mmode[i] = modem[i];
        data->omega_nm[i] = omega_nm[i];
        data->phase_nm[i] = phase_nm[i];
        data->amplitude_nm[i] = amplitude_nm[i];
 
        err = interp1Dcomp_setup(&data->alpha_nm[i], &alpha[i*nrho],
                                 nrho, NATURALBC, rhomin, rhomax);
        if(err) {
            print_err("Error: Failed to initialize splines.\n");
            return err;
        }
        err = interp1Dcomp_setup(&data->phi_nm[i], &phi[i*nrho],
                                 nrho, NATURALBC, rhomin, rhomax);
        if(err) {
            print_err("Error: Failed to initialize splines.\n");
            return err;
        }
    }
 
    /* Print some sanity check on data */
    print_out(VERBOSE_IO, "\nMHD (stationary) input\n");
    print_out(VERBOSE_IO, "Grid: nrho = %4.d rhomin = %3.3f rhomax = %3.3f\n",
              nrho, data->rho_min, data->rho_max);
 
    print_out(VERBOSE_IO, "\nModes:\n");
    for(int i = 0; i < nmode; i++) {
        print_out(VERBOSE_IO,
                  "(n,m) = (%2.d,%2.d) Amplitude = %3.3g Frequency = %3.3g"
                  " Phase = %3.3g\n",
                  data->nmode[i], data->mmode[i], data->amplitude_nm[i],
                  data->omega_nm[i], data->phase_nm[i]);
    }
 
    return err;
}

void mhd_stat_free(mhd_stat_data* data) {
    for(int i = 0; i < data->n_modes; i++) {
        free(data->phi_nm[i].c);
        free(data->alpha_nm[i].c);
    }
    free(data->phi_nm);
    free(data->alpha_nm);
    free(data->nmode);
    free(data->mmode);
    free(data->phase_nm);
    free(data->omega_nm);
    free(data->amplitude_nm);
}

void mhd_stat_offload(mhd_stat_data* data) {
    //TODO: Implement
}

a5err mhd_stat_eval(real mhd_dmhd[10], real r, real phi, real z, real t,
                    int includemode, boozer_data* boozerdata,
                    mhd_stat_data* mhddata, B_field_data* Bdata) {
 
    a5err err = 0;
 
    real ptz[12];
    int isinside;
    if(!err) {
        err = boozer_eval_psithetazeta(ptz, &isinside, r, phi, z, Bdata,
                                       boozerdata);
    }
    real rho[2];
    if(!err && isinside) {
        err = B_field_eval_rho(rho, ptz[0], Bdata);
    }
 
    /* Initialize values */
    for(int i=0; i<10; i++) {
        mhd_dmhd[i] = 0;
    }
 
    int interperr = 0;
    for(int i = 0; i < mhddata->n_modes; i++){
        if( includemode != MHD_INCLUDE_ALL && includemode != i ) { continue; }
        /* Get interpolated values */
        real a_da[3], phi_dphi[3];
        interperr += interp1Dcomp_eval_df(a_da, &(mhddata->alpha_nm[i]),
                                          rho[0]);
        interperr += interp1Dcomp_eval_df(phi_dphi, &(mhddata->phi_nm[i]),
                                          rho[0]);
 
        /* The interpolation returns dx/drho but we require dx/dpsi.
         * The second order derivatives are not needed anywhere */
        a_da[1]     *= rho[1];
        phi_dphi[1] *= rho[1];
 
        /* These are used frequently, so store them in separate variables */
        real mhdarg = mhddata->nmode[i] * ptz[8]
            - mhddata->mmode[i] * ptz[4]
            - mhddata->omega_nm[i] * t
            + mhddata->phase_nm[i];
        real sinmhd = sin(mhdarg);
        real cosmhd = cos(mhdarg);
 
        /* Sum over modes to get alpha, phi */
        mhd_dmhd[0] +=     a_da[0] * mhddata->amplitude_nm[i] * cosmhd;
        mhd_dmhd[5] += phi_dphi[0] * mhddata->amplitude_nm[i] * cosmhd;
 
        /* Time derivatives */
        mhd_dmhd[1] +=     a_da[0] * mhddata->amplitude_nm[i]
            * mhddata->omega_nm[i] * sinmhd;
        mhd_dmhd[6] += phi_dphi[0] * mhddata->amplitude_nm[i]
            * mhddata->omega_nm[i] * sinmhd;
 
        /* R component of gradients */
        mhd_dmhd[2] += mhddata->amplitude_nm[i]
            * (  a_da[1] * ptz[1] * cosmhd
               + a_da[0] * mhddata->mmode[i] * ptz[5] * sinmhd
               - a_da[0] * mhddata->nmode[i] * ptz[9] * sinmhd);
        mhd_dmhd[7] += mhddata->amplitude_nm[i]
            * (   phi_dphi[1] * ptz[1] * cosmhd
                + phi_dphi[0] * mhddata->mmode[i] * ptz[5] * sinmhd
                - phi_dphi[0] * mhddata->nmode[i] * ptz[9] * sinmhd);
 
        /* phi component of gradients */
        mhd_dmhd[3] += (1/r) * mhddata->amplitude_nm[i]
            * (  a_da[1] * ptz[2] * cosmhd
               + a_da[0] * mhddata->mmode[i] * ptz[6]  * sinmhd
               - a_da[0] * mhddata->nmode[i] * ptz[10] * sinmhd);
        mhd_dmhd[8] += (1/r) * mhddata->amplitude_nm[i]
            * (   phi_dphi[1] * ptz[2] * cosmhd
                + phi_dphi[0] * mhddata->mmode[i] * ptz[6]  * sinmhd
                - phi_dphi[0] * mhddata->nmode[i] * ptz[10] * sinmhd);
 
        /* z component of gradients */
        mhd_dmhd[4] += mhddata->amplitude_nm[i]
            * (   a_da[1] * ptz[3] * cosmhd
                + a_da[0] * mhddata->mmode[i] * ptz[7]  * sinmhd
                - a_da[0] * mhddata->nmode[i] * ptz[11] * sinmhd);
        mhd_dmhd[9] += mhddata->amplitude_nm[i]
            * (   phi_dphi[1] * ptz[3] * cosmhd
                + phi_dphi[0] * mhddata->mmode[i] * ptz[7]  * sinmhd
                - phi_dphi[0] * mhddata->nmode[i] * ptz[11] * sinmhd);
    }
 
    /* Omit evaluation if point outside the boozer or mhd grid. */
    if(!isinside || interperr) {
        for(int i=0; i<10; i++) {
            mhd_dmhd[i] = 0;
        }
    }
    return err;
}

a5err mhd_stat_perturbations(real pert_field[7], real r, real phi, real z,
                             real t, int pertonly, int includemode,
                             boozer_data* boozerdata, mhd_stat_data* mhddata,
                             B_field_data* Bdata) {
    a5err err = 0;
    real mhd_dmhd[10];
    if(!err) {
        err = mhd_stat_eval(mhd_dmhd, r, phi, z, t, includemode, boozerdata,
                            mhddata, Bdata);
    }
    /*  see example of curl evaluation in step_gc_rk4.c, ydot_gc*/
    real B_dB[15];
    if(!err) {
        err = B_field_eval_B_dB(B_dB, r, phi, z, t, Bdata);
    }
 
    if(!err) {
        real B[3];
        B[0] = B_dB[0];
        B[1] = B_dB[4];
        B[2] = B_dB[8];
 
        real curlB[3];
        curlB[0] = B_dB[10]/r - B_dB[7];
        curlB[1] = B_dB[3] - B_dB[9];
        curlB[2] = (B[1] - B_dB[2])/r + B_dB[5];
 
        real gradalpha[3];
        gradalpha[0] = mhd_dmhd[2];
        gradalpha[1] = mhd_dmhd[3];
        gradalpha[2] = mhd_dmhd[4];
 
        real gradalphacrossB[3];
 
        math_cross(gradalpha, B, gradalphacrossB);
 
        pert_field[0] = mhd_dmhd[0]*curlB[0] + gradalphacrossB[0];
        pert_field[1] = mhd_dmhd[0]*curlB[1] + gradalphacrossB[1];
        pert_field[2] = mhd_dmhd[0]*curlB[2] + gradalphacrossB[2];
 
        pert_field[3] = -mhd_dmhd[7] - B[0]*mhd_dmhd[1];
        pert_field[4] = -mhd_dmhd[8] - B[1]*mhd_dmhd[1];
        pert_field[5] = -mhd_dmhd[9] - B[2]*mhd_dmhd[1];
        pert_field[6] = mhd_dmhd[5];
 
        if(!pertonly) {
            pert_field[0] += B[0];
            pert_field[1] += B[1];
            pert_field[2] += B[2];
        }
    }
 
    return err;
}