plasma_1Dt.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_1Dt.h"
 
 
int plasma_1Dt_init_offload(plasma_1Dt_offload_data* offload_data,
                            real** offload_array) {
 
    int n_rho = offload_data->n_rho;
    int n_time = offload_data->n_time;
    int n_ions = offload_data->n_species -1;
    print_out(VERBOSE_IO, "\n1D plasma profiles (P_1Dt)\n");
    print_out(VERBOSE_IO,
              "Min rho = %1.2le, Max rho = %1.2le,"
              " Number of rho grid points = %d\n",
              (*offload_array)[0], (*offload_array)[n_rho-1], n_rho);
    print_out(VERBOSE_IO,
              "Min time = %1.2le, Max time = %1.2le,"
              " Number of time points = %d\n",
              (*offload_array)[n_rho], (*offload_array)[n_rho+n_time-1],n_time);
    print_out(VERBOSE_IO,
              "Number of ion species = %d\n", n_ions);
    print_out(VERBOSE_IO,
              "Species Z/A  charge [e]/mass [amu] "
              "Density [m^-3] at Min/Max rho(t=t0)"
              "  Temperature [eV] at Min/Max rho(t=t0)\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)[2*n_rho*n_time + n_time + n_rho*(3+i)],
                  (*offload_array)[2*n_rho*n_time + n_time + n_rho*(4+i) - 1],
                  (*offload_array)[n_time + n_rho*2] / CONST_E,
                  (*offload_array)[n_time + 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)[2*n_rho*n_time + n_time + n_rho*1],
              (*offload_array)[2*n_rho*n_time + n_time + n_rho*2 - 1],
              (*offload_array)[n_time + n_rho*1] / CONST_E,
              (*offload_array)[n_time + n_rho*2 - 1] / CONST_E);
    real quasineutrality = 0;
    for(int k = 0; k <n_rho; k++) {
        real ele_qdens =
            (*offload_array)[2*n_rho*n_time + n_time + n_rho + k] * CONST_E;
        real ion_qdens = 0;
        for(int i=0; i < n_ions; i++) {
            int idx = 2*n_rho*n_time + n_time + n_rho * (2+1) + k;
            ion_qdens +=
                (*offload_array)[idx] * 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_1Dt_free_offload(plasma_1Dt_offload_data* offload_data,
                             real** offload_array) {
    free(*offload_array);
    *offload_array = NULL;
}
 
void plasma_1Dt_init(plasma_1Dt_data* pls_data,
                     plasma_1Dt_offload_data* offload_data,
                     real* offload_array) {
 
    pls_data->n_rho = offload_data->n_rho;
    pls_data->n_time = offload_data->n_time;
    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->rho  = &offload_array[0];
    pls_data->time = &offload_array[pls_data->n_rho];
    pls_data->temp = &offload_array[pls_data->n_rho+pls_data->n_time];
    pls_data->dens = &offload_array[pls_data->n_rho+pls_data->n_time
                                    +2*pls_data->n_rho*pls_data->n_time];
}
 
a5err plasma_1Dt_eval_temp(real* temp, real rho, real t, int species,
                           plasma_1Dt_data* pls_data) {
 
    real temp_dens[MAX_SPECIES], temp_temp[MAX_SPECIES];
 
    a5err err = plasma_1Dt_eval_densandtemp(temp_dens, temp_temp, rho, t,
                                            pls_data);
 
    *temp = temp_temp[species];
 
    return err;
 
}
 
a5err plasma_1Dt_eval_dens(real* dens, real rho, real t, int species,
                           plasma_1Dt_data* pls_data) {
    real temp_dens[MAX_SPECIES], temp_temp[MAX_SPECIES];
 
    a5err err = plasma_1Dt_eval_densandtemp(temp_dens, temp_temp, rho, t,
                                            pls_data);
 
    *dens = temp_dens[species];
 
    return err;
 
}
 
a5err plasma_1Dt_eval_densandtemp(real* dens, real* temp, real rho, real t,
                                  plasma_1Dt_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]);
 
        int i_time = 0;
        while(i_time < pls_data->n_time-1 && pls_data->time[i_time] <= t) {
            i_time++;
        }
        i_time--;
 
        real t_time = (t - pls_data->time[i_time])
                 / (pls_data->time[i_time+1] - pls_data->time[i_time]);
 
        if(i_time < 0) {
            /* time < t[0], use first profile */
            i_time = 0;
            t_time = 0;
        }
        else if(i_time >= pls_data->n_time-2) {
            /* time > t[n_time-1], use last profile */
            i_time = pls_data->n_time-2;
            t_time = 1;
        }
 
        for(int i = 0; i < pls_data->n_species; i++) {
            real p11, p12, p21, p22, p1, p2;
 
            p11 = pls_data->dens[i_time*pls_data->n_species*pls_data->n_rho
                                 + i*pls_data->n_rho
                                 + i_rho];
            p12 = pls_data->dens[i_time*pls_data->n_species*pls_data->n_rho
                                 + i*pls_data->n_rho
                                 + i_rho + 1];
            p21 = pls_data->dens[(i_time+1)*pls_data->n_species*pls_data->n_rho
                                 + i*pls_data->n_rho
                                 + i_rho];
            p22 = pls_data->dens[(i_time+1)*pls_data->n_species*pls_data->n_rho
                                 + i*pls_data->n_rho
                                 + i_rho + 1];
 
            p1 = p11 + t_rho * (p12 - p11);
            p2 = p21 + t_rho * (p22 - p21);
 
            dens[i] = p1 + t_time * (p2 - p1);
 
            if(i < 2) {
                /* Electron and ion temperature */
                p11 = pls_data->temp[i_time*2*pls_data->n_rho
                                     + i*pls_data->n_rho
                                     +i_rho];
                p12 = pls_data->temp[i_time*2*pls_data->n_rho
                                     + i*pls_data->n_rho
                                     + i_rho + 1];
                p21 = pls_data->temp[(i_time+1)*2*pls_data->n_rho
                                     + i*pls_data->n_rho
                                     + i_rho];
                p22 = pls_data->temp[(i_time+1)*2*pls_data->n_rho
                                     + i*pls_data->n_rho
                                     + i_rho + 1];
 
                p1 = p11 + t_rho * (p12 - p11);
                p2 = p21 + t_rho * (p22 - p21);
 
                temp[i] = p1 + t_time * (p2 - p1);
            }
            else {
                /* Temperature is same for all ion species */
                temp[i] = temp[1];
            }
        }
    }
 
    return err;
}