B_3DS.c

Go to the documentation of this file.

 
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "../math.h"
#include "../ascot5.h"
#include "../error.h"
#include "../print.h"
#include "B_3DS.h"
#include "../spline/interp.h"
 
int psigrid_n_r;     
    int psigrid_n_z;     
    real psigrid_r_min;  
    real psigrid_r_max;  
    real psigrid_z_min;  
    real psigrid_z_max;  
    int Bgrid_n_r;       
    int Bgrid_n_z;       
    real Bgrid_r_min;    
    real Bgrid_r_max;    
    real Bgrid_z_min;    
    real Bgrid_z_max;    
    int Bgrid_n_phi;     
    real Bgrid_phi_min;  
    real Bgrid_phi_max;  
    real psi0;           
    real psi1;           
    real axis_r;         
    real axis_z;         
int B_3DS_init(B_3DS_data* data,
               int p_n_r, real p_r_min, real p_r_max,
               int p_n_z, real p_z_min, real p_z_max,
               int b_n_r, real b_r_min, real b_r_max,
               int b_n_phi, real b_phi_min, real b_phi_max,
               int b_n_z, real b_z_min, real b_z_max,
               real axis_r, real axis_z, real psi0, real psi1,
               real* psi, real* B_r, real* B_phi, real* B_z) {
 
    int err = 0;
    data->psi0   = psi0;
    data->psi1   = psi1;
    data->axis_r = axis_r;
    data->axis_z = axis_z;
 
    /* Set up the splines */
    err = interp2Dcomp_setup(&data->psi, psi, p_n_r, p_n_z,
                             NATURALBC, NATURALBC,
                             p_r_min, p_r_max, p_z_min, p_z_max);
    if(err) {
        print_err("Error: Failed to initialize splines.\n");
        return 1;
    }
    err = interp3Dcomp_setup(&data->B_r, B_r, b_n_r, b_n_phi, b_n_z,
                             NATURALBC, PERIODICBC, NATURALBC,
                             b_r_min, b_r_max, b_phi_min, b_phi_max,
                             b_z_min, b_z_max);
    if(err) {
        print_err("Error: Failed to initialize splines.\n");
        return 1;
    }
    err = interp3Dcomp_setup(&data->B_phi, B_phi, b_n_r, b_n_phi, b_n_z,
                             NATURALBC, PERIODICBC, NATURALBC,
                             b_r_min, b_r_max, b_phi_min, b_phi_max,
                             b_z_min, b_z_max);
    if(err) {
        print_err("Error: Failed to initialize splines.\n");
        return 1;
    }
    err = interp3Dcomp_setup(&data->B_z, B_z, b_n_r, b_n_phi, b_n_z,
                             NATURALBC, PERIODICBC, NATURALBC,
                             b_r_min, b_r_max, b_phi_min, b_phi_max,
                             b_z_min, b_z_max);
    if(err) {
        print_err("Error: Failed to initialize splines.\n");
        return 1;
    }
 
    /* Evaluate psi and magnetic field on axis for checks */
    real psival[1], Bval[3];
    err = B_3DS_eval_psi(psival, data->axis_r, 0, data->axis_z, data);
    err = B_3DS_eval_B(Bval, data->axis_r, 0, data->axis_z, data);
    if(err) {
        print_err("Error: Initialization failed.\n");
        return err;
    }
 
    /* Print some sanity check on data */
    printf("\n3D magnetic field (B_3DS)\n");
    print_out(VERBOSE_IO, "Psi-grid: nR = %4.d Rmin = %3.3f m Rmax = %3.3f m\n",
              p_n_r, p_r_min, p_r_max);
    print_out(VERBOSE_IO, "      nz = %4.d zmin = %3.3f m zmax = %3.3f m\n",
              p_n_z, p_z_min, p_z_max);
    print_out(VERBOSE_IO, "B-grid: nR = %4.d Rmin = %3.3f m Rmax = %3.3f m\n",
              b_n_r, b_r_min, b_r_max);
    print_out(VERBOSE_IO, "      nz = %4.d zmin = %3.3f m zmax = %3.3f m\n",
              b_n_z, b_z_min, b_z_max);
    print_out(VERBOSE_IO, "nphi = %4.d phimin = %3.3f deg phimax = %3.3f deg\n",
              b_n_phi, math_rad2deg(b_phi_min), math_rad2deg(b_phi_max));
    print_out(VERBOSE_IO, "Psi at magnetic axis (%1.3f m, %1.3f m)\n",
              data->axis_r, data->axis_z);
    print_out(VERBOSE_IO, "%3.3f (evaluated)\n%3.3f (given)\n",
              psival[0], data->psi0);
    print_out(VERBOSE_IO, "Magnetic field on axis:\n"
              "B_R = %3.3f B_phi = %3.3f B_z = %3.3f\n",
              Bval[0], Bval[1], Bval[2]);
 
    return err;
}
 
void B_3DS_free(B_3DS_data* data) {
    free(data->psi.c);
    free(data->B_r.c);
    free(data->B_phi.c);
    free(data->B_z.c);
}
 
void B_3DS_offload(B_3DS_data* data) {
    GPU_MAP_TO_DEVICE(
        data->psi, data->B_r, data->B_phi, data->B_z, \
        data->psi.c[0:data->psi.n_x*data->psi.n_y*NSIZE_COMP2D], \
        data->B_r.c[0:data->B_r.n_x*data->B_r.n_y*data->B_r.n_z*NSIZE_COMP3D], \
        data->B_phi.c[0:data->B_phi.n_x*data->B_phi.n_y*data->B_phi.n_z*NSIZE_COMP3D], \
        data->B_z.c[0:data->B_z.n_x*data->B_z.n_y*data->B_z.n_z*NSIZE_COMP3D]
    )
}
 
a5err B_3DS_eval_psi(real* psi, real r, real phi, real z,
                   B_3DS_data* Bdata) {
    a5err err = 0;
    int interperr = 0; /* If error happened during interpolation */
    interperr += interp2Dcomp_eval_f(&psi[0], &Bdata->psi, r, z);
 
    if(interperr) {
        err = error_raise( ERR_INPUT_EVALUATION, __LINE__, EF_B_3DS );
    }
 
    return err;
}
 
a5err B_3DS_eval_psi_dpsi(real psi_dpsi[4], real r, real phi, real z,
                   B_3DS_data* Bdata) {
    a5err err = 0;
    int interperr = 0;
    real psi_dpsi_temp[6];
 
    interperr += interp2Dcomp_eval_df(psi_dpsi_temp, &Bdata->psi, r, z);
 
    psi_dpsi[0] = psi_dpsi_temp[0];
    psi_dpsi[1] = psi_dpsi_temp[1];
    psi_dpsi[2] = 0;
    psi_dpsi[3] = psi_dpsi_temp[2];
 
    if(interperr) {
        err = error_raise( ERR_INPUT_EVALUATION, __LINE__, EF_B_3DS );
    }
 
    return err;
}
 
a5err B_3DS_eval_rho_drho(real rho_drho[4], real r, real phi, real z,
                          B_3DS_data* Bdata) {
    int interperr = 0; /* If error happened during interpolation */
    real psi_dpsi[6];
 
    interperr += interp2Dcomp_eval_df(psi_dpsi, &Bdata->psi, r, z);
 
    if(interperr) {
        return error_raise( ERR_INPUT_EVALUATION, __LINE__, EF_B_3DS );
    }
 
    /* Check that the values seem valid */
    real delta = Bdata->psi1 - Bdata->psi0;
    if( (psi_dpsi[0] - Bdata->psi0) / delta < 0 ) {
         return error_raise( ERR_INPUT_UNPHYSICAL, __LINE__, EF_B_3DS );
    }
 
    /* Normalize psi to get rho */
    rho_drho[0] = sqrt(fabs((psi_dpsi[0] - Bdata->psi0) / delta));
 
    rho_drho[1] = psi_dpsi[1] / (2*delta*rho_drho[0]);
    rho_drho[2] = 0;
    rho_drho[3] = psi_dpsi[2] / (2*delta*rho_drho[0]);
 
    return 0;
}
 
a5err B_3DS_eval_B(real B[3], real r, real phi, real z, B_3DS_data* Bdata) {
    a5err err = 0;
    int interperr = 0;
 
    interperr += interp3Dcomp_eval_f(&B[0], &Bdata->B_r, r, phi, z);
    interperr += interp3Dcomp_eval_f(&B[1], &Bdata->B_phi, r, phi, z);
    interperr += interp3Dcomp_eval_f(&B[2], &Bdata->B_z, r, phi, z);
 
    /* Test for B field interpolation error */
    if(interperr) {
        err = error_raise( ERR_INPUT_EVALUATION, __LINE__, EF_B_3DS );
    }
 
    if(!err) {
        real psi_dpsi[6];
        interperr += interp2Dcomp_eval_df(psi_dpsi, &Bdata->psi, r, z);
 
        B[0] = B[0] - psi_dpsi[2]/r;
        B[2] = B[2] + psi_dpsi[1]/r;
 
        /* Test for psi interpolation error */
        if(interperr) {
            err = error_raise( ERR_INPUT_EVALUATION, __LINE__, EF_B_3DS );
        }
    }
 
    /* Check that magnetic field seems valid */
    int check = 0;
    check += ((B[0]*B[0] + B[1]*B[1] + B[2]*B[2]) == 0);
    if(!err && check) {
        err = error_raise( ERR_INPUT_UNPHYSICAL, __LINE__, EF_B_3DS );
    }
 
    return err;
}
 
a5err B_3DS_eval_B_dB(real B_dB[12], real r, real phi, real z,
                      B_3DS_data* Bdata) {
    a5err err = 0;
    int interperr = 0; /* If error happened during interpolation */
    real B_dB_temp[10];
 
    interperr += interp3Dcomp_eval_df(B_dB_temp, &Bdata->B_r, r, phi, z);
    B_dB[0] = B_dB_temp[0];
    B_dB[1] = B_dB_temp[1];
    B_dB[2] = B_dB_temp[2];
    B_dB[3] = B_dB_temp[3];
 
    interperr += interp3Dcomp_eval_df(B_dB_temp, &Bdata->B_phi, r, phi, z);
    B_dB[4] = B_dB_temp[0];
    B_dB[5] = B_dB_temp[1];
    B_dB[6] = B_dB_temp[2];
    B_dB[7] = B_dB_temp[3];
 
    interperr += interp3Dcomp_eval_df(B_dB_temp, &Bdata->B_z, r, phi, z);
    B_dB[8] = B_dB_temp[0];
    B_dB[9] = B_dB_temp[1];
    B_dB[10] = B_dB_temp[2];
    B_dB[11] = B_dB_temp[3];
 
    /* Test for B field interpolation error */
    if(interperr) {
        err = error_raise( ERR_INPUT_EVALUATION, __LINE__, EF_B_3DS );
    }
 
    if(!err) {
        real psi_dpsi[6];
        interperr += interp2Dcomp_eval_df(psi_dpsi, &Bdata->psi, r, z);
 
        B_dB[0] = B_dB[0] - psi_dpsi[2]/r;
        B_dB[1] = B_dB[1] + psi_dpsi[2]/(r*r)-psi_dpsi[5]/r;
        B_dB[3] = B_dB[3] - psi_dpsi[4]/r;
        B_dB[8] = B_dB[8] + psi_dpsi[1]/r;
        B_dB[9] = B_dB[9] - psi_dpsi[1]/(r*r) + psi_dpsi[3]/r;
        B_dB[11] = B_dB[11] + psi_dpsi[5]/r;
 
        /* Test for psi interpolation error */
        if(interperr) {
            err = error_raise( ERR_INPUT_EVALUATION, __LINE__, EF_B_3DS );
        }
    }
 
    /* Check that magnetic field seems valid */
    int check = 0;
    check += ((B_dB[0]*B_dB[0] + B_dB[4]*B_dB[4] + B_dB[8]*B_dB[8]) == 0);
    if(!err && check) {
        err = error_raise( ERR_INPUT_UNPHYSICAL, __LINE__, EF_B_3DS );
    }
 
    return err;
}
 
a5err B_3DS_get_axis_rz(real rz[2], B_3DS_data* Bdata) {
    a5err err = 0;
    rz[0] = Bdata->axis_r;
    rz[1] = Bdata->axis_z;
    return err;
}