plasma_1DS.c

Go to the documentation of this file.

 
#include <stdlib.h>
#include <math.h>
#include <float.h>
#include "../ascot5.h"
#include "../print.h"
#include "../error.h"
#include "../consts.h"
#include "../spline/interp.h"
#include "plasma_1DS.h"
 
#define PLASMA_1DS_NONEG 1
 
#define PLASMA_1DS_LOG 1
 
#define PLASMA_1DS_SQRT 2
 
int plasma_1DS_init_offload(plasma_1DS_offload_data* offload_data,
                            real** offload_array) {
 
    /* Spline initialization */
    int err = 0;
    int n_rho     = offload_data->n_rho;
    int n_species = offload_data->n_species;
 
    /* Allocate enough space for two temperature and n_species density arrays */
    real* coeff_array = (real*) malloc((2 + n_species) * NSIZE_COMP1D
                                       * n_rho * sizeof(real));
 
    for(int i=0; i< ((2 + n_species)*n_rho); i++) {
#if PLASMA_1DS_NONEG == PLASMA_1DS_LOG
        (*offload_array)[i] = log( (*offload_array)[i]);
#elif PLASMA_1DS_NONEG == PLASMA_1DS_SQRT
        (*offload_array)[i] = sqrt( (*offload_array)[i] );
#endif
    }
 
    /* Evaluate spline coefficients */
 
    /* Te */
    err += interp1Dcomp_init_coeff(
        coeff_array + 0*n_rho*NSIZE_COMP1D,
        *offload_array + 0*n_rho,
        offload_data->n_rho, NATURALBC,
        offload_data->rho_min, offload_data->rho_max);
 
    /* Ti */
    err += interp1Dcomp_init_coeff(
        coeff_array + 1*n_rho*NSIZE_COMP1D,
        *offload_array + 1*n_rho,
        offload_data->n_rho, NATURALBC,
        offload_data->rho_min, offload_data->rho_max);
 
    /* Densities */
    for(int i=0; i < n_species; i++) {
        err += interp1Dcomp_init_coeff(
            coeff_array + (2 +i)*n_rho*NSIZE_COMP1D,
            *offload_array + (2 + i)*n_rho,
            offload_data->n_rho, NATURALBC,
            offload_data->rho_min, offload_data->rho_max);
    }
 
    if(err) {
        free(coeff_array);
        return err;
    }
 
    free(*offload_array);
    *offload_array = coeff_array;
    offload_data->offload_array_length = (2 + n_species) * NSIZE_COMP1D * n_rho;
 
    /* Initialize the data so that we can perform sanity checks */
    real temp0, temp1, dens0, dens1;
    plasma_1DS_data plsdata;
    plasma_1DS_init(&plsdata, offload_data, *offload_array);
 
    int n_ions = n_species - 1;
    print_out(VERBOSE_IO, "\n1D plasma profiles (P_1DS)\n");
    print_out(VERBOSE_IO,
              "Min rho = %1.2le, Max rho = %1.2le,"
              " Number of rho grid points = %d,"
              " Number of ion species = %d\n",
              offload_data->rho_min, offload_data->rho_max, n_rho, n_ions);
    print_out(VERBOSE_IO,
              "Species Z/A  charge [e]/mass [amu] Density [m^-3] at Min/Max rho"
              "    Temperature [eV] at Min/Max rho\n");
    for(int i=0; i < n_ions; i++) {
        plasma_1DS_eval_temp(&temp0, offload_data->rho_min, i+1, &plsdata);
        plasma_1DS_eval_temp(&temp1, offload_data->rho_max, i+1, &plsdata);
        plasma_1DS_eval_dens(&dens0, offload_data->rho_min, i+1, &plsdata);
        plasma_1DS_eval_dens(&dens1, offload_data->rho_max, i+1, &plsdata);
        print_out(VERBOSE_IO,
                  " %3d  /%3d   %3d  /%7.3f             %1.2le/%1.2le     "
                  "           %1.2le/%1.2le       \n",
                  offload_data->znum[i], offload_data->anum[i],
                  (int)round(offload_data->charge[i+1]/CONST_E),
                  offload_data->mass[i+1]/CONST_U,
                  dens0, dens1, temp0 / CONST_E, temp1 / CONST_E);
    }
 
    plasma_1DS_eval_temp(&temp0, offload_data->rho_min, 0, &plsdata);
    plasma_1DS_eval_temp(&temp1, offload_data->rho_max, 0, &plsdata);
    plasma_1DS_eval_dens(&dens0, offload_data->rho_min, 0, &plsdata);
    plasma_1DS_eval_dens(&dens1, offload_data->rho_max, 0, &plsdata);
    print_out(VERBOSE_IO,
              "[electrons]  %3d  /%7.3f             %1.2le/%1.2le          "
              "      %1.2le/%1.2le       \n", -1, CONST_M_E/CONST_U,
              dens0, dens1, temp0 / CONST_E, temp1 / CONST_E);
    real quasineutrality = 0;
    for(int k = 0; k <n_rho; k++) {
        real rho = offload_data->rho_min
            + k * (offload_data->rho_max - offload_data->rho_min) / (n_rho - 1);
 
        real ele_qdens;
        plasma_1DS_eval_dens(&ele_qdens, rho, 0, &plsdata);
        real ion_qdens = 0;
        for(int i=0; i < n_ions; i++) {
            plasma_1DS_eval_dens(&dens0, rho, i+1, &plsdata);
            ion_qdens += dens0;
        }
        quasineutrality = fmax( quasineutrality,
                                fabs( 1 - ion_qdens / ele_qdens ) );
    }
    print_out(VERBOSE_IO, "Quasi-neutrality is (electron / ion charge density)"
              " %.2f\n", 1+quasineutrality);
 
    return 0;
}
 
void plasma_1DS_free_offload(plasma_1DS_offload_data* offload_data,
                             real** offload_array) {
    free(*offload_array);
    *offload_array = NULL;
}
 
void plasma_1DS_init(plasma_1DS_data* plasma_data,
                     plasma_1DS_offload_data* offload_data,
                     real* offload_array) {
    plasma_data->n_species = offload_data->n_species;
 
    for(int i = 0; i < plasma_data->n_species; i++) {
        plasma_data->mass[i]   = offload_data->mass[i];
        plasma_data->charge[i] = offload_data->charge[i];
        plasma_data->znum[i]   = offload_data->znum[i];
        plasma_data->anum[i]   = offload_data->anum[i];
    }
 
    int n_rho = offload_data->n_rho;
    interp1Dcomp_init_spline(&(plasma_data->temp[0]),
                             &(offload_array[0*n_rho]),
                             offload_data->n_rho, NATURALBC,
                             offload_data->rho_min,
                             offload_data->rho_max);
    interp1Dcomp_init_spline(&(plasma_data->temp[1]),
                             &(offload_array[1*n_rho*NSIZE_COMP1D]),
                             offload_data->n_rho, NATURALBC,
                             offload_data->rho_min,
                             offload_data->rho_max);
 
    for(int i=0; i<offload_data->n_species; i++) {
        interp1Dcomp_init_spline(&(plasma_data->dens[i]),
                                 &(offload_array[(2+i)*n_rho*NSIZE_COMP1D]),
                                 offload_data->n_rho, NATURALBC,
                                 offload_data->rho_min,
                                 offload_data->rho_max);
    }
}
 
a5err plasma_1DS_eval_temp(real* temp, real rho, int species,
                           plasma_1DS_data* plasma_data) {
    int interperr = 0;
    interperr += interp1Dcomp_eval_f(temp, &plasma_data->temp[species>0], rho);
 
    a5err err = 0;
    if(interperr) {
        err = error_raise( ERR_INPUT_EVALUATION, __LINE__, EF_PLASMA_1DS );
    }
 
#if PLASMA_1DS_NONEG == PLASMA_1DS_LOG
    *temp = exp(*temp);
#elif PLASMA_1DS_NONEG == PLASMA_1DS_SQRT
    *temp = (*temp) * (*temp);
#endif
    if(!err && *temp < 0){
        err = error_raise( ERR_INPUT_EVALUATION, __LINE__, EF_PLASMA_1DS );
    }
    return err;
}
 
a5err plasma_1DS_eval_dens(real* dens, real rho, int species,
                           plasma_1DS_data* plasma_data) {
 
    int interperr = 0;
    interperr += interp1Dcomp_eval_f(dens, &plasma_data->dens[species], rho);
 
    a5err err = 0;
    if(interperr) {
        err = error_raise( ERR_INPUT_EVALUATION, __LINE__, EF_PLASMA_1DS );
    }
 
#if PLASMA_1DS_NONEG == PLASMA_1DS_LOG
    *dens = exp(*dens);
#elif PLASMA_1DS_NONEG == PLASMA_1DS_SQRT
    *dens = (*dens) * (*dens);
#endif
    if(!err && *dens < 0){
        err = error_raise( ERR_INPUT_EVALUATION, __LINE__, EF_PLASMA_1DS );
    }
    return err;
}
 
a5err plasma_1DS_eval_densandtemp(real* dens, real* temp, real rho,
                                  plasma_1DS_data* plasma_data) {
    int interperr = 0;
 
    /* Evaluate electron temperature and density */
    interperr += interp1Dcomp_eval_f(&temp[0], &plasma_data->temp[0], rho);
    interperr += interp1Dcomp_eval_f(&dens[0], &plasma_data->dens[0], rho);
 
    /* Evaluate ion temperature (same for all ions) and densities */
    interperr += interp1Dcomp_eval_f(&temp[1], &plasma_data->temp[1], rho);
    for(int i=1; i<plasma_data->n_species; i++) {
        temp[i]    = temp[1];
        interperr += interp1Dcomp_eval_f(&dens[i], &plasma_data->dens[i], rho);
    }
 
    a5err err = 0;
    if(interperr) {
        err = error_raise( ERR_INPUT_EVALUATION, __LINE__, EF_PLASMA_1DS );
    }
 
#if PLASMA_1DS_NONEG == PLASMA_1DS_LOG
    for(int i=0; i<plasma_data->n_species; i++) {
        dens[i] = exp(dens[i]);
        temp[i] = exp(temp[i]);
    }
#elif PLASMA_1DS_NONEG == PLASMA_1DS_SQRT
    for(int i=0; i<plasma_data->n_species; i++) {
        dens[i] = dens[i]*dens[i];
        temp[i] = temp[i]*temp[i];
    }
#endif
    for(int i=0; i<plasma_data->n_species; i++) {
        if(!err && (dens[i] < 0 || temp[i] < 0) ) {
            err = error_raise( ERR_INPUT_EVALUATION, __LINE__,
                                EF_PLASMA_1DS );
        }
    }
    return err;
}