plasma_1D.c

Go to the documentation of this file.

 
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "../ascot5.h"
#include "../error.h"
#include "../consts.h"
#include "../print.h"
#include "plasma_1D.h"
 
 
int plasma_1D_init_offload(plasma_1D_offload_data* offload_data,
                           real** offload_array) {
 
    int n_rho = offload_data->n_rho;
    int n_ions = offload_data->n_species -1;
    print_out(VERBOSE_IO, "\n1D plasma profiles (P_1D)\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_array)[0], (*offload_array)[n_rho-1], 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++) {
        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,
                  (*offload_array)[n_rho*(4+i)],
                  (*offload_array)[n_rho*(5+i) - 1],
                  (*offload_array)[n_rho*2] / CONST_E,
                  (*offload_array)[n_rho*3-1] / CONST_E);
    }
    print_out(VERBOSE_IO,
              "[electrons]  %3d  /%7.3f             %1.2le/%1.2le          "
              "      %1.2le/%1.2le       \n",
              -1, CONST_M_E/CONST_U,
              (*offload_array)[n_rho*3],
              (*offload_array)[n_rho*4 - 1],
              (*offload_array)[n_rho] / CONST_E,
              (*offload_array)[n_rho*2-1] / CONST_E);
    real quasineutrality = 0;
    for(int k = 0; k <n_rho; k++) {
        real ele_qdens = (*offload_array)[n_rho*3 + k] * CONST_E;
        real ion_qdens = 0;
        for(int i=0; i < n_ions; i++) {
            ion_qdens +=
                (*offload_array)[n_rho*(4+i) + k] * offload_data->charge[i+1];
        }
        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_1D_free_offload(plasma_1D_offload_data* offload_data,
                            real** offload_array) {
    free(*offload_array);
    *offload_array = NULL;
}
 
void plasma_1D_init(plasma_1D_data* pls_data,
                    plasma_1D_offload_data* offload_data,
                    real* offload_array) {
 
    pls_data->n_rho = offload_data->n_rho;
    pls_data->n_species = offload_data->n_species;
 
    for(int i = 0; i < pls_data->n_species; i++) {
        pls_data->mass[i]   = offload_data->mass[i];
        pls_data->charge[i] = offload_data->charge[i];
        pls_data->znum[i]   = offload_data->znum[i];
        pls_data->anum[i]   = offload_data->anum[i];
    }
    pls_data->rho  = &offload_array[0];
    pls_data->temp = &offload_array[pls_data->n_rho];
    pls_data->dens = &offload_array[pls_data->n_rho*3];
}
 
a5err plasma_1D_eval_temp(real* temp, real rho, int species,
                          plasma_1D_data* pls_data) {
 
    a5err err = 0;
    if(rho < pls_data->rho[0]) {
        err = error_raise( ERR_INPUT_EVALUATION, __LINE__, EF_PLASMA_1D );
    }
    else if(rho >= pls_data->rho[pls_data->n_rho-1]) {
        err = error_raise( ERR_INPUT_EVALUATION, __LINE__, EF_PLASMA_1D );
    }
    else {
        int i_rho = 0;
        while(i_rho < pls_data->n_rho - 1 && pls_data->rho[i_rho] <= rho) {
            i_rho++;
        }
        i_rho--;
        real t_rho = (rho - pls_data->rho[i_rho])
            / (pls_data->rho[i_rho+1] - pls_data->rho[i_rho]);
 
        real p1 = pls_data->temp[species*pls_data->n_rho + i_rho];
        real p2 = pls_data->temp[species*pls_data->n_rho + i_rho+1];
        temp[0] = p1 + t_rho * (p2 - p1);
    }
 
    return err;
}
 
a5err plasma_1D_eval_dens(real* dens, real rho, int species,
                          plasma_1D_data* pls_data) {
 
    a5err err = 0;
    if(rho < pls_data->rho[0]) {
        err = error_raise( ERR_INPUT_EVALUATION, __LINE__, EF_PLASMA_1D );
    }
    else if(rho >= pls_data->rho[pls_data->n_rho-1]) {
        err = error_raise( ERR_INPUT_EVALUATION, __LINE__, EF_PLASMA_1D );
    }
    else {
        int i_rho = 0;
        while(i_rho < pls_data->n_rho - 1 && pls_data->rho[i_rho] <= rho) {
            i_rho++;
        }
        i_rho--;
        real t_rho = (rho - pls_data->rho[i_rho])
                 / (pls_data->rho[i_rho+1] - pls_data->rho[i_rho]);
 
        real p1 = pls_data->dens[species*pls_data->n_rho + i_rho];
        real p2 = pls_data->dens[species*pls_data->n_rho + i_rho+1];
        dens[0] = p1 + t_rho * (p2 - p1);
    }
 
    return err;
}
 
a5err plasma_1D_eval_densandtemp(real* dens, real* temp, real rho,
                                 plasma_1D_data* pls_data) {
    a5err err = 0;
    if(rho < pls_data->rho[0]) {
        err = error_raise( ERR_INPUT_EVALUATION, __LINE__, EF_PLASMA_1D );
    }
    else if(rho >= pls_data->rho[pls_data->n_rho-1]) {
        err = error_raise( ERR_INPUT_EVALUATION, __LINE__, EF_PLASMA_1D );
    }
    else {
        int i_rho = 0;
        while(i_rho < pls_data->n_rho-1 && pls_data->rho[i_rho] <= rho) {
            i_rho++;
        }
        i_rho--;
 
        real t_rho = (rho - pls_data->rho[i_rho])
                 / (pls_data->rho[i_rho+1] - pls_data->rho[i_rho]);
 
        real p1, p2;
        for(int i = 0; i < pls_data->n_species; i++) {
            p1 = pls_data->dens[i*pls_data->n_rho + i_rho];
            p2 = pls_data->dens[i*pls_data->n_rho + i_rho+1];
            dens[i] = p1 + t_rho * (p2 - p1);
 
            if(i < 2) {
                /* Electron and ion temperature */
                p1 = pls_data->temp[i*pls_data->n_rho + i_rho];
                p2 = pls_data->temp[i*pls_data->n_rho + i_rho+1];
                temp[i] = p1 + t_rho * (p2 - p1);
            }
            else {
                /* Temperature is same for all ion species */
                temp[i] = temp[1];
            }
        }
    }
 
    return err;
}