B_STS.c

Go to the documentation of this file.

 
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include "../math.h"
#include "../ascot5.h"
#include "../error.h"
#include "../print.h"
#include "../consts.h"
#include "B_STS.h"
#include "../linint/linint.h"
#include "../spline/interp.h"
 
int B_STS_init(B_STS_data* data,
               int p_n_r, real p_r_min, real p_r_max,
               int p_n_phi, real p_phi_min, real p_phi_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,
               int naxis, real axis_min, real axis_max,
               real* axis_r, real* axis_z, real psi0, real psi1,
               real* psi, real* B_r, real* B_phi, real* B_z) {
 
    /* Spline initialization. */
    int err = 0;
    data->psi0 = psi0;
    data->psi1 = psi1;
    err = interp3Dcomp_setup(&data->psi, psi, p_n_r, p_n_phi, p_n_z,
                             NATURALBC, PERIODICBC, NATURALBC,
                             p_r_min, p_r_max, p_phi_min, p_phi_max,
                             p_z_min, p_z_max);
    if(err) {
        print_err("Error: Failed to initialize splines.\n");
        return 1;
    }
    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;
    }
    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;
    }
    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;
    }
 
    real* c1 = (real*)malloc(naxis*sizeof(real));
    real* c2 = (real*)malloc(naxis*sizeof(real));
    for(int i = 0; i < naxis; i++) {
        c1[i] = axis_r[i];
        c2[i] = axis_z[i];
    }
    linint1D_init(&data->axis_r, c1, naxis, PERIODICBC, axis_min, axis_max);
    linint1D_init(&data->axis_z, c2, naxis, PERIODICBC, axis_min, axis_max);
 
    /* Evaluate psi and magnetic field on axis for checks */
    real psival[1], Bval[3], axis[2];
    err += B_STS_get_axis_rz(axis, data, 0);
    err += B_STS_eval_psi(psival, axis[0], 0, axis[1], data);
    err += B_STS_eval_B(Bval, axis[0], 0, axis[1], data);
    if(err) {
        print_err("Error: Initialization failed.\n");
        return err;
    }
 
    printf("\nStellarator magnetic field (B_STS)\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, "nphi = %4.d phimin = %3.3f deg phimax = %3.3f deg\n",
              p_n_phi, math_rad2deg(p_phi_min), math_rad2deg(p_phi_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 (phi=0) (%1.3f m, %1.3f m)\n",
              axis[0], axis[1]);
    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 0;
}
 
void B_STS_free(B_STS_data* data) {
    free(data->psi.c);
    free(data->B_r.c);
    free(data->B_phi.c);
    free(data->B_z.c);
}
 
void B_STS_offload(B_STS_data* data) {
    GPU_MAP_TO_DEVICE(
        data->axis_r, data->axis_r.c[0:data->axis_r.n_x], \
        data->axis_z, data->axis_z.c[0:data->axis_z.n_x], \
        data->psi, data->B_r, data->B_phi, data->B_z, \
        data->psi.c[0:data->psi.n_x*data->psi.n_y*data->psi.n_z*NSIZE_COMP3D], \
        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_STS_eval_psi(real* psi, real r, real phi, real z,
                     B_STS_data* Bdata) {
    a5err err = 0;
    int interperr = 0; /* If error happened during interpolation */
 
    interperr += interp3Dcomp_eval_f(&psi[0], &Bdata->psi, r, phi, z);
 
#ifdef B_STS_CLAMP_RHO_NONNEGATIVE
    if ( psi[0] < Bdata->psi0 ){
        psi[0] = Bdata->psi0;
    }
#endif
 
    /* Test for psi interpolation error */
    if(interperr) {
        return error_raise( ERR_INPUT_EVALUATION, __LINE__, EF_B_STS );
    }
 
    return err;
}
 
a5err B_STS_eval_psi_dpsi(real psi_dpsi[4], real r, real phi, real z,
                          B_STS_data* Bdata) {
    a5err err = 0;
    int interperr = 0; /* If error happened during interpolation */
    real psi_dpsi_temp[10];
 
    interperr += interp3Dcomp_eval_df(psi_dpsi_temp, &Bdata->psi, r, phi, z);
 
 
    psi_dpsi[0] = psi_dpsi_temp[0];
    psi_dpsi[1] = psi_dpsi_temp[1];
    psi_dpsi[2] = psi_dpsi_temp[2];
    psi_dpsi[3] = psi_dpsi_temp[3];
 
#ifdef B_STS_CLAMP_RHO_NONNEGATIVE
    if ( psi_dpsi_temp[0] < Bdata->psi0 ){
        psi_dpsi[0] = Bdata->psi0;
        psi_dpsi[1] = 0.0;
        psi_dpsi[2] = 0.0;
        psi_dpsi[3] = 0.0;
    }
#endif
 
 
    if(interperr) {
        return error_raise( ERR_INPUT_EVALUATION, __LINE__, EF_B_STS );
    }
 
    return err;
}
 
a5err B_STS_eval_rho_drho(real rho_drho[4], real r, real phi, real z,
                          B_STS_data* Bdata) {
    a5err err = 0;
    real psi_dpsi[4];
 
    err = B_STS_eval_psi_dpsi(psi_dpsi, r, phi, z, Bdata);
    if(err){
        return error_raise( ERR_INPUT_UNPHYSICAL, __LINE__, EF_B_STS );
    }
 
    /* Check that the values seem valid */
    real delta = Bdata->psi1 - Bdata->psi0;
    if( (psi_dpsi[0] - Bdata->psi0) / delta <= 0 ) {
#ifdef B_STS_CLAMP_RHO_NONNEGATIVE
        /* Set rho and the partial derivatives to zero, because otherwise one
           would need to divide by rho, which is zero. Of course, this problem
           persists when B_STS_CLAMP_RHO_NONNEGATIVE is not defined and
           the evaluation happens exactly at rho=0.0 */
        rho_drho[0] = 0.0;
        rho_drho[1] = 0.0;
        rho_drho[2] = 0.0;
        rho_drho[3] = 0.0;
        return err;
#else
        return error_raise( ERR_INPUT_UNPHYSICAL, __LINE__, EF_B_STS );
#endif
    }
 
    /* Normalize psi to get rho */
    rho_drho[0] = sqrt( (psi_dpsi[0] - Bdata->psi0) / delta );
 
    rho_drho[1] = psi_dpsi[1] / (2*delta*rho_drho[0]);
    rho_drho[2] = psi_dpsi[2] / (2*delta*rho_drho[0]);
    rho_drho[3] = psi_dpsi[3] / (2*delta*rho_drho[0]);
 
    return err;
}
 
a5err B_STS_eval_B(real B[3], real r, real phi, real z, B_STS_data* Bdata) {
    a5err err = 0;
    int interperr = 0; /* If error happened during interpolation */
 
    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) {
        return error_raise( ERR_INPUT_EVALUATION, __LINE__, EF_B_STS );
    }
 
    /* 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) {
        return error_raise( ERR_INPUT_UNPHYSICAL, __LINE__, EF_B_STS );
    }
 
    return err;
 
}
 
a5err B_STS_eval_B_dB(real B_dB[12], real r, real phi, real z,
                      B_STS_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) {
        return error_raise( ERR_INPUT_EVALUATION, __LINE__, EF_B_STS );
    }
 
    /* 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) {
        return error_raise( ERR_INPUT_UNPHYSICAL, __LINE__, EF_B_STS );
    }
 
    return 0;
}
 
a5err B_STS_get_axis_rz(real rz[2], B_STS_data* Bdata, real phi) {
    a5err err = 0;
 
    int interperr = 0; /* If error happened during interpolation */
    interperr += linint1D_eval_f(&rz[0], &Bdata->axis_r, phi);
    interperr += linint1D_eval_f(&rz[1], &Bdata->axis_z, phi);
    if(interperr) {
        err = error_raise( ERR_INPUT_EVALUATION, __LINE__, EF_B_STS );
    }
    return err;
}