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_offload(B_STS_offload_data* offload_data, real** offload_array) {
 
    /* Spline initialization. */
    int err = 0;
    int psi_size = offload_data->psigrid_n_r * offload_data->psigrid_n_z
                   * offload_data->psigrid_n_phi;
    int B_size   = offload_data->Bgrid_n_r   * offload_data->Bgrid_n_z
                   * offload_data->Bgrid_n_phi;
    int axis_size = offload_data->n_axis;
 
    /* Allocate enough space to store four 3D arrays and axis data */
    real* coeff_array = (real*) malloc( (3*NSIZE_COMP3D*B_size
                                         + NSIZE_COMP3D*psi_size
                                         + 2*axis_size)*sizeof(real));
    real* B_r    = &(coeff_array[0*B_size*NSIZE_COMP3D]);
    real* B_phi  = &(coeff_array[1*B_size*NSIZE_COMP3D]);
    real* B_z    = &(coeff_array[2*B_size*NSIZE_COMP3D]);
    real* psi    = &(coeff_array[3*B_size*NSIZE_COMP3D]);
    real* axis_r = &(coeff_array[(3*B_size + psi_size)*NSIZE_COMP3D]);
    real* axis_z = &(coeff_array[(3*B_size + psi_size)*NSIZE_COMP3D
                                 + axis_size]);
 
    err += interp3Dcomp_init_coeff(
        psi, *offload_array + 3*B_size,
        offload_data->psigrid_n_r, offload_data->psigrid_n_phi,
        offload_data->psigrid_n_z,
        NATURALBC, PERIODICBC, NATURALBC,
        offload_data->psigrid_r_min,   offload_data->psigrid_r_max,
        offload_data->psigrid_phi_min, offload_data->psigrid_phi_max,
        offload_data->psigrid_z_min,   offload_data->psigrid_z_max);
 
    err += interp3Dcomp_init_coeff(
        B_r, *offload_array + 0*B_size,
        offload_data->Bgrid_n_r, offload_data->Bgrid_n_phi,
        offload_data->Bgrid_n_z,
        NATURALBC, PERIODICBC, NATURALBC,
        offload_data->Bgrid_r_min,   offload_data->Bgrid_r_max,
        offload_data->Bgrid_phi_min, offload_data->Bgrid_phi_max,
        offload_data->Bgrid_z_min,   offload_data->Bgrid_z_max);
 
    err += interp3Dcomp_init_coeff(
        B_phi, *offload_array + 1*B_size,
        offload_data->Bgrid_n_r, offload_data->Bgrid_n_phi,
        offload_data->Bgrid_n_z,
        NATURALBC, PERIODICBC, NATURALBC,
        offload_data->Bgrid_r_min,   offload_data->Bgrid_r_max,
        offload_data->Bgrid_phi_min, offload_data->Bgrid_phi_max,
        offload_data->Bgrid_z_min,   offload_data->Bgrid_z_max);
 
    err += interp3Dcomp_init_coeff(
        B_z, *offload_array + 2*B_size,
        offload_data->Bgrid_n_r, offload_data->Bgrid_n_phi,
        offload_data->Bgrid_n_z,
        NATURALBC, PERIODICBC, NATURALBC,
        offload_data->Bgrid_r_min,   offload_data->Bgrid_r_max,
        offload_data->Bgrid_phi_min, offload_data->Bgrid_phi_max,
        offload_data->Bgrid_z_min,   offload_data->Bgrid_z_max);
 
    if(err) {
        print_err("Error: Failed to initialize splines.\n");
        return err;
    }
 
    for(int i = 0; i < axis_size; i++) {
        axis_r[i] = (*offload_array)[3*B_size + psi_size + i];
        axis_z[i] = (*offload_array)[3*B_size + psi_size + axis_size + i];
    }
 
    /* Re-allocate the offload array and store spline coefficients there */
    free(*offload_array);
    *offload_array = coeff_array;
    offload_data->offload_array_length = 3*NSIZE_COMP3D*B_size
                                         + NSIZE_COMP3D*psi_size
                                         + 2*axis_size;
 
    /* Evaluate psi and magnetic field on axis for checks */
    B_STS_data Bdata;
    B_STS_init(&Bdata, offload_data, *offload_array);
    real psival[1], Bval[3], axis[2];
    err += B_STS_get_axis_rz(axis, &Bdata, 0);
    err += B_STS_eval_psi(psival, axis[0], 0, axis[1],
                         &Bdata);
    err += B_STS_eval_B(Bval, axis[0], 0, axis[1],
                       &Bdata);
    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",
              offload_data->psigrid_n_r,
              offload_data->psigrid_r_min, offload_data->psigrid_r_max);
    print_out(VERBOSE_IO, "      nz = %4.d zmin = %3.3f m zmax = %3.3f m\n",
              offload_data->psigrid_n_z,
              offload_data->psigrid_z_min, offload_data->psigrid_z_max);
    print_out(VERBOSE_IO, "nphi = %4.d phimin = %3.3f deg phimax = %3.3f deg\n",
              offload_data->psigrid_n_phi,
              math_rad2deg(offload_data->psigrid_phi_min),
              math_rad2deg(offload_data->psigrid_phi_max));
    print_out(VERBOSE_IO, "B-grid: nR = %4.d Rmin = %3.3f m Rmax = %3.3f m\n",
              offload_data->Bgrid_n_r,
              offload_data->Bgrid_r_min, offload_data->Bgrid_r_max);
    print_out(VERBOSE_IO, "      nz = %4.d zmin = %3.3f m zmax = %3.3f m\n",
              offload_data->Bgrid_n_z,
              offload_data->Bgrid_z_min, offload_data->Bgrid_z_max);
    print_out(VERBOSE_IO, "nphi = %4.d phimin = %3.3f deg phimax = %3.3f deg\n",
              offload_data->Bgrid_n_phi,
              math_rad2deg(offload_data->Bgrid_phi_min),
              math_rad2deg(offload_data->Bgrid_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], offload_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_offload(B_STS_offload_data* offload_data, real** offload_array) {
    free(*offload_array);
    *offload_array = NULL;
}
 
void B_STS_init(B_STS_data* Bdata, B_STS_offload_data* offload_data,
                real* offload_array) {
 
    int B_size = offload_data->Bgrid_n_r * offload_data->Bgrid_n_z
        * offload_data->Bgrid_n_phi*NSIZE_COMP3D;
    int psi_size = offload_data->psigrid_n_r * offload_data->psigrid_n_z
        * offload_data->psigrid_n_phi*NSIZE_COMP3D;
    int axis_size = offload_data->n_axis;
 
    /* Initialize target data struct */
    Bdata->psi0 = offload_data->psi0;
    Bdata->psi1 = offload_data->psi1;
 
 
    /* Initialize spline structs from the coefficients */
    interp3Dcomp_init_spline(&Bdata->B_r, &(offload_array[0*B_size]),
                             offload_data->Bgrid_n_r,
                             offload_data->Bgrid_n_phi,
                             offload_data->Bgrid_n_z,
                             NATURALBC, PERIODICBC, NATURALBC,
                             offload_data->Bgrid_r_min,
                             offload_data->Bgrid_r_max,
                             offload_data->Bgrid_phi_min,
                             offload_data->Bgrid_phi_max,
                             offload_data->Bgrid_z_min,
                             offload_data->Bgrid_z_max);
 
    interp3Dcomp_init_spline(&Bdata->B_phi, &(offload_array[1*B_size]),
                             offload_data->Bgrid_n_r,
                             offload_data->Bgrid_n_phi,
                             offload_data->Bgrid_n_z,
                             NATURALBC, PERIODICBC, NATURALBC,
                             offload_data->Bgrid_r_min,
                             offload_data->Bgrid_r_max,
                             offload_data->Bgrid_phi_min,
                             offload_data->Bgrid_phi_max,
                             offload_data->Bgrid_z_min,
                             offload_data->Bgrid_z_max);
 
    interp3Dcomp_init_spline(&Bdata->B_z, &(offload_array[2*B_size]),
                             offload_data->Bgrid_n_r,
                             offload_data->Bgrid_n_phi,
                             offload_data->Bgrid_n_z,
                             NATURALBC, PERIODICBC, NATURALBC,
                             offload_data->Bgrid_r_min,
                             offload_data->Bgrid_r_max,
                             offload_data->Bgrid_phi_min,
                             offload_data->Bgrid_phi_max,
                             offload_data->Bgrid_z_min,
                             offload_data->Bgrid_z_max);
 
    interp3Dcomp_init_spline(&Bdata->psi, &(offload_array[3*B_size]),
                             offload_data->psigrid_n_r,
                             offload_data->psigrid_n_phi,
                             offload_data->psigrid_n_z,
                             NATURALBC, PERIODICBC, NATURALBC,
                             offload_data->psigrid_r_min,
                             offload_data->psigrid_r_max,
                             offload_data->psigrid_phi_min,
                             offload_data->psigrid_phi_max,
                             offload_data->psigrid_z_min,
                             offload_data->psigrid_z_max);
 
    linint1D_init(&Bdata->axis_r,
                  &(offload_array[3*B_size + psi_size]),
                  offload_data->n_axis, PERIODICBC,
                  offload_data->axis_min, offload_data->axis_max);
 
    linint1D_init(&Bdata->axis_z,
                  &(offload_array[3*B_size + psi_size + axis_size]),
                  offload_data->n_axis, PERIODICBC,
                  offload_data->axis_min, offload_data->axis_max);
}
 
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;
}