ASCOT5
Loading...
Searching...
No Matches
B_3DS.c
Go to the documentation of this file.
1
34#include <stdlib.h>
35#include <stdio.h>
36#include <math.h>
37#include "../math.h"
38#include "../ascot5.h"
39#include "../error.h"
40#include "../print.h"
41#include "B_3DS.h"
42#include "../spline/interp.h"
43
93int B_3DS_init_offload(B_3DS_offload_data* offload_data, real** offload_array) {
94
95 /* Spline initialization. */
96 int err = 0;
97 int psi_size = offload_data->psigrid_n_r * offload_data->psigrid_n_z;
98 int B_size = offload_data->Bgrid_n_r * offload_data->Bgrid_n_z
99 * offload_data->Bgrid_n_phi;
100
101 /* Allocate enough space to store three 3D arrays and one 2D array */
102 real* coeff_array = (real*) malloc( (3*NSIZE_COMP3D*B_size
103 + NSIZE_COMP2D*psi_size)*sizeof(real));
104 real* B_r = &(coeff_array[0*B_size*NSIZE_COMP3D]);
105 real* B_phi = &(coeff_array[1*B_size*NSIZE_COMP3D]);
106 real* B_z = &(coeff_array[2*B_size*NSIZE_COMP3D]);
107 real* psi = &(coeff_array[3*B_size*NSIZE_COMP3D]);
108
109 err += interp2Dcomp_init_coeff(
110 psi, *offload_array + 3*B_size,
111 offload_data->psigrid_n_r, offload_data->psigrid_n_z,
112 NATURALBC, NATURALBC,
113 offload_data->psigrid_r_min, offload_data->psigrid_r_max,
114 offload_data->psigrid_z_min, offload_data->psigrid_z_max);
115
116 err += interp3Dcomp_init_coeff(
117 B_r, *offload_array + 0*B_size,
118 offload_data->Bgrid_n_r, offload_data->Bgrid_n_phi,
119 offload_data->Bgrid_n_z,
120 NATURALBC, PERIODICBC, NATURALBC,
121 offload_data->Bgrid_r_min, offload_data->Bgrid_r_max,
122 offload_data->Bgrid_phi_min, offload_data->Bgrid_phi_max,
123 offload_data->Bgrid_z_min, offload_data->Bgrid_z_max);
124
125 err += interp3Dcomp_init_coeff(
126 B_phi, *offload_array + 1*B_size,
127 offload_data->Bgrid_n_r, offload_data->Bgrid_n_phi,
128 offload_data->Bgrid_n_z,
129 NATURALBC, PERIODICBC, NATURALBC,
130 offload_data->Bgrid_r_min, offload_data->Bgrid_r_max,
131 offload_data->Bgrid_phi_min, offload_data->Bgrid_phi_max,
132 offload_data->Bgrid_z_min, offload_data->Bgrid_z_max);
133
134 err += interp3Dcomp_init_coeff(
135 B_z, *offload_array + 2*B_size,
136 offload_data->Bgrid_n_r, offload_data->Bgrid_n_phi,
137 offload_data->Bgrid_n_z,
138 NATURALBC, PERIODICBC, NATURALBC,
139 offload_data->Bgrid_r_min, offload_data->Bgrid_r_max,
140 offload_data->Bgrid_phi_min, offload_data->Bgrid_phi_max,
141 offload_data->Bgrid_z_min, offload_data->Bgrid_z_max);
142
143 if(err) {
144 print_err("Error: Failed to initialize splines.\n");
145 return err;
146 }
147
148 /* Re-allocate the offload array and store spline coefficients there */
149 free(*offload_array);
150 *offload_array = coeff_array;
151 offload_data->offload_array_length = NSIZE_COMP2D*psi_size
152 + NSIZE_COMP3D*B_size*3;
153
154 /* Evaluate psi and magnetic field on axis for checks */
155 B_3DS_data Bdata;
156 B_3DS_init(&Bdata, offload_data, *offload_array);
157 real psival[1], Bval[3];
158 err = B_3DS_eval_psi(psival, offload_data->axis_r, 0, offload_data->axis_z,
159 &Bdata);
160 err = B_3DS_eval_B(Bval, offload_data->axis_r, 0, offload_data->axis_z,
161 &Bdata);
162 if(err) {
163 print_err("Error: Initialization failed.\n");
164 return err;
165 }
166
167 /* Print some sanity check on data */
168 printf("\n3D magnetic field (B_3DS)\n");
169 print_out(VERBOSE_IO, "Psi-grid: nR = %4.d Rmin = %3.3f m Rmax = %3.3f m\n",
170 offload_data->psigrid_n_r,
171 offload_data->psigrid_r_min, offload_data->psigrid_r_max);
172 print_out(VERBOSE_IO, " nz = %4.d zmin = %3.3f m zmax = %3.3f m\n",
173 offload_data->psigrid_n_z,
174 offload_data->psigrid_z_min, offload_data->psigrid_z_max);
175 print_out(VERBOSE_IO, "B-grid: nR = %4.d Rmin = %3.3f m Rmax = %3.3f m\n",
176 offload_data->Bgrid_n_r,
177 offload_data->Bgrid_r_min, offload_data->Bgrid_r_max);
178 print_out(VERBOSE_IO, " nz = %4.d zmin = %3.3f m zmax = %3.3f m\n",
179 offload_data->Bgrid_n_z,
180 offload_data->Bgrid_z_min, offload_data->Bgrid_z_max);
181 print_out(VERBOSE_IO, "nphi = %4.d phimin = %3.3f deg phimax = %3.3f deg\n",
182 offload_data->Bgrid_n_phi,
183 math_rad2deg(offload_data->Bgrid_phi_min),
184 math_rad2deg(offload_data->Bgrid_phi_max));
185 print_out(VERBOSE_IO, "Psi at magnetic axis (%1.3f m, %1.3f m)\n",
186 offload_data->axis_r, offload_data->axis_z);
187 print_out(VERBOSE_IO, "%3.3f (evaluated)\n%3.3f (given)\n",
188 psival[0], offload_data->psi0);
189 print_out(VERBOSE_IO, "Magnetic field on axis:\n"
190 "B_R = %3.3f B_phi = %3.3f B_z = %3.3f\n",
191 Bval[0], Bval[1], Bval[2]);
192
193 return err;
194}
195
202void B_3DS_free_offload(B_3DS_offload_data* offload_data,
203 real** offload_array) {
204 free(*offload_array);
205 *offload_array = NULL;
206}
207
215void B_3DS_init(B_3DS_data* Bdata, B_3DS_offload_data* offload_data,
216 real* offload_array) {
217
218 int B_size = NSIZE_COMP3D * offload_data->Bgrid_n_r
219 * offload_data->Bgrid_n_z * offload_data->Bgrid_n_phi;
220
221 /* Initialize target data struct */
222 Bdata->psi0 = offload_data->psi0;
223 Bdata->psi1 = offload_data->psi1;
224 Bdata->axis_r = offload_data->axis_r;
225 Bdata->axis_z = offload_data->axis_z;
226
227 /* Initialize spline structs from the coefficients */
228 interp3Dcomp_init_spline(&Bdata->B_r, &(offload_array[0*B_size]),
229 offload_data->Bgrid_n_r,
230 offload_data->Bgrid_n_phi,
231 offload_data->Bgrid_n_z,
232 NATURALBC, PERIODICBC, NATURALBC,
233 offload_data->Bgrid_r_min,
234 offload_data->Bgrid_r_max,
235 offload_data->Bgrid_phi_min,
236 offload_data->Bgrid_phi_max,
237 offload_data->Bgrid_z_min,
238 offload_data->Bgrid_z_max);
239
240 interp3Dcomp_init_spline(&Bdata->B_phi, &(offload_array[1*B_size]),
241 offload_data->Bgrid_n_r,
242 offload_data->Bgrid_n_phi,
243 offload_data->Bgrid_n_z,
244 NATURALBC, PERIODICBC, NATURALBC,
245 offload_data->Bgrid_r_min,
246 offload_data->Bgrid_r_max,
247 offload_data->Bgrid_phi_min,
248 offload_data->Bgrid_phi_max,
249 offload_data->Bgrid_z_min,
250 offload_data->Bgrid_z_max);
251
252 interp3Dcomp_init_spline(&Bdata->B_z, &(offload_array[2*B_size]),
253 offload_data->Bgrid_n_r,
254 offload_data->Bgrid_n_phi,
255 offload_data->Bgrid_n_z,
256 NATURALBC, PERIODICBC, NATURALBC,
257 offload_data->Bgrid_r_min,
258 offload_data->Bgrid_r_max,
259 offload_data->Bgrid_phi_min,
260 offload_data->Bgrid_phi_max,
261 offload_data->Bgrid_z_min,
262 offload_data->Bgrid_z_max);
263
264 interp2Dcomp_init_spline(&Bdata->psi, &(offload_array[3*B_size]),
265 offload_data->psigrid_n_r,
266 offload_data->psigrid_n_z,
267 NATURALBC, NATURALBC,
268 offload_data->psigrid_r_min,
269 offload_data->psigrid_r_max,
270 offload_data->psigrid_z_min,
271 offload_data->psigrid_z_max);
272}
273
285a5err B_3DS_eval_psi(real* psi, real r, real phi, real z,
286 B_3DS_data* Bdata) {
287 a5err err = 0;
288 int interperr = 0; /* If error happened during interpolation */
289 interperr += interp2Dcomp_eval_f(&psi[0], &Bdata->psi, r, z);
290
291 if(interperr) {
292 err = error_raise( ERR_INPUT_EVALUATION, __LINE__, EF_B_3DS );
293 }
294
295 return err;
296}
297
309a5err B_3DS_eval_psi_dpsi(real psi_dpsi[4], real r, real phi, real z,
310 B_3DS_data* Bdata) {
311 a5err err = 0;
312 int interperr = 0;
313 real psi_dpsi_temp[6];
314
315 interperr += interp2Dcomp_eval_df(psi_dpsi_temp, &Bdata->psi, r, z);
316
317 psi_dpsi[0] = psi_dpsi_temp[0];
318 psi_dpsi[1] = psi_dpsi_temp[1];
319 psi_dpsi[2] = 0;
320 psi_dpsi[3] = psi_dpsi_temp[2];
321
322 if(interperr) {
323 err = error_raise( ERR_INPUT_EVALUATION, __LINE__, EF_B_3DS );
324 }
325
326 return err;
327}
328
340a5err B_3DS_eval_rho_drho(real rho_drho[4], real r, real phi, real z,
341 B_3DS_data* Bdata) {
342 int interperr = 0; /* If error happened during interpolation */
343 real psi_dpsi[6];
344
345 interperr += interp2Dcomp_eval_df(psi_dpsi, &Bdata->psi, r, z);
346
347 if(interperr) {
348 return error_raise( ERR_INPUT_EVALUATION, __LINE__, EF_B_3DS );
349 }
350
351 /* Check that the values seem valid */
352 real delta = Bdata->psi1 - Bdata->psi0;
353 if( (psi_dpsi[0] - Bdata->psi0) / delta < 0 ) {
354 return error_raise( ERR_INPUT_UNPHYSICAL, __LINE__, EF_B_3DS );
355 }
356
357 /* Normalize psi to get rho */
358 rho_drho[0] = sqrt(fabs((psi_dpsi[0] - Bdata->psi0) / delta));
359
360 rho_drho[1] = psi_dpsi[1] / (2*delta*rho_drho[0]);
361 rho_drho[2] = 0;
362 rho_drho[3] = psi_dpsi[2] / (2*delta*rho_drho[0]);
363
364 return 0;
365}
366
378a5err B_3DS_eval_B(real B[3], real r, real phi, real z, B_3DS_data* Bdata) {
379 a5err err = 0;
380 int interperr = 0;
381
382 interperr += interp3Dcomp_eval_f(&B[0], &Bdata->B_r, r, phi, z);
383 interperr += interp3Dcomp_eval_f(&B[1], &Bdata->B_phi, r, phi, z);
384 interperr += interp3Dcomp_eval_f(&B[2], &Bdata->B_z, r, phi, z);
385
386 /* Test for B field interpolation error */
387 if(interperr) {
388 err = error_raise( ERR_INPUT_EVALUATION, __LINE__, EF_B_3DS );
389 }
390
391 if(!err) {
392 real psi_dpsi[6];
393 interperr += interp2Dcomp_eval_df(psi_dpsi, &Bdata->psi, r, z);
394
395 B[0] = B[0] - psi_dpsi[2]/r;
396 B[2] = B[2] + psi_dpsi[1]/r;
397
398 /* Test for psi interpolation error */
399 if(interperr) {
400 err = error_raise( ERR_INPUT_EVALUATION, __LINE__, EF_B_3DS );
401 }
402 }
403
404 /* Check that magnetic field seems valid */
405 int check = 0;
406 check += ((B[0]*B[0] + B[1]*B[1] + B[2]*B[2]) == 0);
407 if(!err && check) {
408 err = error_raise( ERR_INPUT_UNPHYSICAL, __LINE__, EF_B_3DS );
409 }
410
411 return err;
412}
413
425a5err B_3DS_eval_B_dB(real B_dB[12], real r, real phi, real z,
426 B_3DS_data* Bdata) {
427 a5err err = 0;
428 int interperr = 0; /* If error happened during interpolation */
429 real B_dB_temp[10];
430
431 interperr += interp3Dcomp_eval_df(B_dB_temp, &Bdata->B_r, r, phi, z);
432 B_dB[0] = B_dB_temp[0];
433 B_dB[1] = B_dB_temp[1];
434 B_dB[2] = B_dB_temp[2];
435 B_dB[3] = B_dB_temp[3];
436
437 interperr += interp3Dcomp_eval_df(B_dB_temp, &Bdata->B_phi, r, phi, z);
438 B_dB[4] = B_dB_temp[0];
439 B_dB[5] = B_dB_temp[1];
440 B_dB[6] = B_dB_temp[2];
441 B_dB[7] = B_dB_temp[3];
442
443 interperr += interp3Dcomp_eval_df(B_dB_temp, &Bdata->B_z, r, phi, z);
444 B_dB[8] = B_dB_temp[0];
445 B_dB[9] = B_dB_temp[1];
446 B_dB[10] = B_dB_temp[2];
447 B_dB[11] = B_dB_temp[3];
448
449 /* Test for B field interpolation error */
450 if(interperr) {
451 err = error_raise( ERR_INPUT_EVALUATION, __LINE__, EF_B_3DS );
452 }
453
454 if(!err) {
455 real psi_dpsi[6];
456 interperr += interp2Dcomp_eval_df(psi_dpsi, &Bdata->psi, r, z);
457
458 B_dB[0] = B_dB[0] - psi_dpsi[2]/r;
459 B_dB[1] = B_dB[1] + psi_dpsi[2]/(r*r)-psi_dpsi[5]/r;
460 B_dB[3] = B_dB[3] - psi_dpsi[4]/r;
461 B_dB[8] = B_dB[8] + psi_dpsi[1]/r;
462 B_dB[9] = B_dB[9] - psi_dpsi[1]/(r*r) + psi_dpsi[3]/r;
463 B_dB[11] = B_dB[11] + psi_dpsi[5]/r;
464
465 /* Test for psi interpolation error */
466 if(interperr) {
467 err = error_raise( ERR_INPUT_EVALUATION, __LINE__, EF_B_3DS );
468 }
469 }
470
471 /* Check that magnetic field seems valid */
472 int check = 0;
473 check += ((B_dB[0]*B_dB[0] + B_dB[4]*B_dB[4] + B_dB[8]*B_dB[8]) == 0);
474 if(!err && check) {
475 err = error_raise( ERR_INPUT_UNPHYSICAL, __LINE__, EF_B_3DS );
476 }
477
478 return err;
479}
480
489a5err B_3DS_get_axis_rz(real rz[2], B_3DS_data* Bdata) {
490 a5err err = 0;
491 rz[0] = Bdata->axis_r;
492 rz[1] = Bdata->axis_z;
493 return err;
494}
B_3DS_free_offload
void B_3DS_free_offload(B_3DS_offload_data *offload_data, real **offload_array)
Free offload array.
Definition B_3DS.c:202
B_3DS_eval_rho_drho
a5err B_3DS_eval_rho_drho(real rho_drho[4], real r, real phi, real z, B_3DS_data *Bdata)
Evaluate normalized poloidal flux rho and its derivatives.
Definition B_3DS.c:340
B_3DS_eval_B
a5err B_3DS_eval_B(real B[3], real r, real phi, real z, B_3DS_data *Bdata)
Evaluate magnetic field.
Definition B_3DS.c:378
B_3DS_init_offload
int B_3DS_init_offload(B_3DS_offload_data *offload_data, real **offload_array)
Initialize magnetic field offload data.
Definition B_3DS.c:93
B_3DS_init
void B_3DS_init(B_3DS_data *Bdata, B_3DS_offload_data *offload_data, real *offload_array)
Initialize magnetic field data struct on target.
Definition B_3DS.c:215
B_3DS_eval_psi_dpsi
a5err B_3DS_eval_psi_dpsi(real psi_dpsi[4], real r, real phi, real z, B_3DS_data *Bdata)
Evaluate poloidal flux psi and its derivatives.
Definition B_3DS.c:309
B_3DS_get_axis_rz
a5err B_3DS_get_axis_rz(real rz[2], B_3DS_data *Bdata)
Return magnetic axis R-coordinate.
Definition B_3DS.c:489
B_3DS_eval_psi
a5err B_3DS_eval_psi(real *psi, real r, real phi, real z, B_3DS_data *Bdata)
Evaluate poloidal flux psi.
Definition B_3DS.c:285
B_3DS_eval_B_dB
a5err B_3DS_eval_B_dB(real B_dB[12], real r, real phi, real z, B_3DS_data *Bdata)
Evaluate magnetic field and its derivatives.
Definition B_3DS.c:425
B_3DS.h
Header file for B_3DS.c.
ascot5.h
Main header file for ASCOT5.
real
double real
Definition ascot5.h:85
error.h
Error module for ASCOT5.
a5err
unsigned long int a5err
Simulation error flag.
Definition error.h:17
EF_B_3DS
@ EF_B_3DS
Definition error.h:36
ERR_INPUT_EVALUATION
@ ERR_INPUT_EVALUATION
Definition error.h:63
ERR_INPUT_UNPHYSICAL
@ ERR_INPUT_UNPHYSICAL
Definition error.h:65
error_raise
static DECLARE_TARGET_SIMD a5err error_raise(error_type type, int line, error_file file)
Raise a new error.
Definition error.h:86
interp.h
Spline interpolation library.
interp3Dcomp_eval_f
DECLARE_TARGET_END a5err interp3Dcomp_eval_f(real *f, interp3D_data *str, real x, real y, real z)
Evaluate interpolated value of 3D scalar field.
Definition interp3Dcomp.c:250
interp3Dcomp_init_coeff
int interp3Dcomp_init_coeff(real *c, real *f, int n_x, int n_y, int n_z, int bc_x, int bc_y, int bc_z, real x_min, real x_max, real y_min, real y_max, real z_min, real z_max)
Calculate tricubic spline interpolation coefficients for 3D data.
Definition interp3Dcomp.c:37
NATURALBC
@ NATURALBC
Definition interp.h:37
PERIODICBC
@ PERIODICBC
Definition interp.h:38
interp2Dcomp_eval_df
DECLARE_TARGET_END a5err interp2Dcomp_eval_df(real *f_df, interp2D_data *str, real x, real y)
Evaluate interpolated value and 1st and 2nd derivatives of 2D field.
Definition interp2Dcomp.c:260
interp2Dcomp_init_coeff
int interp2Dcomp_init_coeff(real *c, real *f, int n_x, int n_y, int bc_x, int bc_y, real x_min, real x_max, real y_min, real y_max)
Calculate bicubic spline interpolation coefficients for scalar 2D data.
Definition interp2Dcomp.c:33
interp3Dcomp_init_spline
void interp3Dcomp_init_spline(interp3D_data *str, real *c, int n_x, int n_y, int n_z, int bc_x, int bc_y, int bc_z, real x_min, real x_max, real y_min, real y_max, real z_min, real z_max)
Initialize a tricubic spline.
Definition interp3Dcomp.c:203
interp3Dcomp_eval_df
DECLARE_TARGET_END a5err interp3Dcomp_eval_df(real *f_df, interp3D_data *str, real x, real y, real z)
Evaluate interpolated value of 3D field and 1st and 2nd derivatives.
Definition interp3Dcomp.c:416
interp2Dcomp_init_spline
void interp2Dcomp_init_spline(interp2D_data *str, real *c, int n_x, int n_y, int bc_x, int bc_y, real x_min, real x_max, real y_min, real y_max)
Initialize a bicubic spline.
Definition interp2Dcomp.c:130
interp2Dcomp_eval_f
DECLARE_TARGET_END a5err interp2Dcomp_eval_f(real *f, interp2D_data *str, real x, real y)
Evaluate interpolated value of a 2D field.
Definition interp2Dcomp.c:168
math.h
Header file for math.c.
math_rad2deg
#define math_rad2deg(a)
Convert radians to degrees.
Definition math.h:110
print.h
Macros for printing console output.
print_out
#define print_out(v,...)
Print to standard output.
Definition print.h:31
VERBOSE_IO
@ VERBOSE_IO
Definition print.h:20
print_err
#define print_err(...)
Print to standard error.
Definition print.h:42
B_3DS_data
3D magnetic field parameters on the target
Definition B_3DS.h:43
B_3DS_offload_data
3D magnetic field parameters on the host
Definition B_3DS.h:15
B_3DS_offload_data::psi1
real psi1
Definition B_3DS.h:34
B_3DS_offload_data::psigrid_z_max
real psigrid_z_max
Definition B_3DS.h:21
B_3DS_offload_data::axis_r
real axis_r
Definition B_3DS.h:35
B_3DS_offload_data::Bgrid_n_z
int Bgrid_n_z
Definition B_3DS.h:24
B_3DS_offload_data::Bgrid_phi_max
real Bgrid_phi_max
Definition B_3DS.h:31
B_3DS_offload_data::psigrid_n_r
int psigrid_n_r
Definition B_3DS.h:16
B_3DS_offload_data::Bgrid_n_phi
int Bgrid_n_phi
Definition B_3DS.h:29
B_3DS_offload_data::psigrid_r_min
real psigrid_r_min
Definition B_3DS.h:18
B_3DS_offload_data::Bgrid_z_max
real Bgrid_z_max
Definition B_3DS.h:28
B_3DS_offload_data::psigrid_n_z
int psigrid_n_z
Definition B_3DS.h:17
B_3DS_offload_data::psi0
real psi0
Definition B_3DS.h:33
B_3DS_offload_data::psigrid_z_min
real psigrid_z_min
Definition B_3DS.h:20
B_3DS_offload_data::Bgrid_n_r
int Bgrid_n_r
Definition B_3DS.h:23
B_3DS_offload_data::axis_z
real axis_z
Definition B_3DS.h:36
B_3DS_offload_data::Bgrid_r_min
real Bgrid_r_min
Definition B_3DS.h:25
B_3DS_offload_data::offload_array_length
int offload_array_length
Definition B_3DS.h:37
B_3DS_offload_data::Bgrid_r_max
real Bgrid_r_max
Definition B_3DS.h:26
B_3DS_offload_data::Bgrid_z_min
real Bgrid_z_min
Definition B_3DS.h:27
B_3DS_offload_data::psigrid_r_max
real psigrid_r_max
Definition B_3DS.h:19
B_3DS_offload_data::Bgrid_phi_min
real Bgrid_phi_min
Definition B_3DS.h:30
Generated on Wed Oct 9 2024 17:13:56 for ASCOT5 by doxygen 1.10.0