wall_2d.c

Go to the documentation of this file.

 
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "../ascot5.h"
#include "../print.h"
#include "wall_2d.h"
 
int wall_2d_init(wall_2d_data* data, int nelements, real* r, real* z) {
 
    data->n = nelements;
    data->wall_r = (real*) malloc( nelements * sizeof(real) );
    data->wall_z = (real*) malloc( nelements * sizeof(real) );
    real rmin = r[0], rmax = r[0];
    real zmin = z[0], zmax = z[0];
    for(int i=0; i < nelements; i++) {
        rmin = fmin(rmin, r[i]);
        rmax = fmax(rmax, r[i]);
        zmin = fmin(zmin, z[i]);
        zmax = fmax(zmax, z[i]);
        data->wall_r[i] = r[i];
        data->wall_z[i] = z[i];
    }
 
    print_out(VERBOSE_IO, "\n2D wall model (wall_2D)\n");
    print_out(VERBOSE_IO, "Number of wall elements = %d,"
              " R extend = [%2.2f, %2.2f], z extend = [%2.2f, %2.2f]\n",
              nelements, rmin, rmax, zmin, zmax);
    return 0;
}
 
void wall_2d_free(wall_2d_data* data) {
    free(data->wall_r);
    free(data->wall_z);
}
 
void wall_2d_offload(wall_2d_data* data) {
    GPU_MAP_TO_DEVICE(
        data->wall_r[0:data->n],
        data->wall_z[0:data->n]
    )
}
 
int wall_2d_inside(real r, real z, wall_2d_data* w) {
    int hits = 0;
    for(int i = 0; i < w->n; i++) {
        real wr1, wr2, wz1, wz2;
        if(i == w->n - 1) {
            wz1 = w->wall_z[i] - z;
            wz2 = w->wall_z[0]   - z;
            wr1 = w->wall_r[i] - r;
            wr2 = w->wall_r[0]   - r;
        } else {
            wz1 = w->wall_z[i]   - z;
            wz2 = w->wall_z[i+1] - z;
            wr1 = w->wall_r[i]   - r;
            wr2 = w->wall_r[i+1] - r;
        }
        if(wz1 * wz2 < 0) {
            real ri = wr1 + (wz1*(wr2-wr1)) / (wz1-wz2);
            if(ri > 0) {
                hits++;
            }
        }
    }
    return hits % 2;
}
 
int wall_2d_hit_wall(real r1, real phi1, real z1, real r2, real phi2, real z2,
                     wall_2d_data* w, real* w_coll) {
    int tile = 0;
    if(!wall_2d_inside(r2, z2, w)) {
        tile = wall_2d_find_intersection(r1, z1, r2, z2, w, w_coll);
    }
    return tile;
}
 
int wall_2d_find_intersection(real r1, real z1, real r2, real z2,
                              wall_2d_data* w, real* w_coll) {
    int tile = 0;
    real t0 = 2.0; // Helper variable to pick the closest intersection
    for(int i=0; i<w->n; i++) {
        real r3, z3, r4, z4;
        if(i == w->n-1) {
            r3 = w->wall_r[i];
            z3 = w->wall_z[i];
            r4 = w->wall_r[0];
            z4 = w->wall_z[0];
        } else {
            r3 = w->wall_r[i];
            z3 = w->wall_z[i];
            r4 = w->wall_r[i+1];
            z4 = w->wall_z[i+1];
        }
 
        real div = (r1 - r2) * (z3 - z4) - (z1 - z2) * (r3 - r4);
        real t   = ( (r1 - r3) * (z3 - z4) - (z1 - z3) * (r3 - r4) ) / div;
        real u   = ( (r1 - r3) * (z1 - z2) - (z1 - z3) * (r1 - r2) ) / div;
        if(0 <= t && t <= 1.0 && 0 <= u && u <= 1.0 && t < t0) {
            t0   = t;
            tile = i + 1;
        }
    }
    *w_coll = t0;
    return tile;
}