ASCOT5
Loading...
Searching...
No Matches
B_GS.c
Go to the documentation of this file.
1
67#include <stdlib.h>
68#include <stdio.h>
69#include <math.h>
70#include "../ascot5.h"
71#include "../consts.h"
72#include "../print.h"
73#include "../error.h"
74#include "B_GS.h"
75
91int B_GS_init_offload(B_GS_offload_data* offload_data, real** offload_array) {
92
93 offload_data->offload_array_length = 0;
94
95 /* Evaluate psi and magnetic field on axis for checks */
96 a5err err = 0;
97 B_GS_data Bdata;
98 B_GS_init(&Bdata, offload_data, *offload_array);
99 real psival[1], Bval[3];
100 err = B_GS_eval_psi(psival, offload_data->raxis, 0, offload_data->zaxis,
101 &Bdata);
102 err = B_GS_eval_B(Bval, offload_data->raxis, 0, offload_data->zaxis,
103 &Bdata);
104 if(err) {
105 print_err("Error: Initialization failed.\n");
106 return err;
107 }
108
109 /* Print some sanity check on data */
110 print_out(VERBOSE_IO,
111 "\nAnalytical tokamak magnetic field (B_GS)\n"
112 "Psi at magnetic axis (%1.3f m, %1.3f m)\n"
113 "%3.3f (evaluated)\n%3.3f (given)\n"
114 "Magnetic field on axis:\n"
115 "B_R = %3.3f B_phi = %3.3f B_z = %3.3f\n"
116 "Number of toroidal field coils %d\n",
117 offload_data->raxis, offload_data->zaxis,
118 psival[0], offload_data->psi0,
119 Bval[0], Bval[1], Bval[2],
120 offload_data->Nripple);
121
122 return 0;
123}
124
131void B_GS_free_offload(B_GS_offload_data* offload_data,
132 real** offload_array) {
133 free(*offload_array);
134 *offload_array = NULL;
135}
136
144void B_GS_init(B_GS_data* Bdata, B_GS_offload_data* offload_data,
145 real* offload_array) {
146 Bdata->R0 = offload_data->R0;
147 Bdata->z0 = offload_data->z0;
148 Bdata->raxis = offload_data->raxis;
149 Bdata->zaxis = offload_data->zaxis;
150 Bdata->B_phi0 = offload_data->B_phi0;
151 Bdata->psi0 = offload_data->psi0;
152 Bdata->psi1 = offload_data->psi1;
153 Bdata->psi_mult = offload_data->psi_mult;
154
155 Bdata->Nripple = offload_data->Nripple;
156 Bdata->a0 = offload_data->a0;
157 Bdata->delta0 = offload_data->delta0;
158 Bdata->alpha0 = offload_data->alpha0;
159
160 Bdata->psi_coeff = offload_data->psi_coeff;
161}
162
174a5err B_GS_eval_psi(real* psi, real r, real phi, real z,
175 B_GS_data* Bdata) {
176 /* Normalize the coordinates */
177 z -= Bdata->z0;
178 r /= Bdata->R0;
179 z /= Bdata->R0;
180
181 real* C = Bdata->psi_coeff;
182
183 /* Precalculated terms used in the components */
184 double r2 = r*r;
185 double r3 = r2*r;
186 double r4 = r3*r;
187 double r5 = r4*r;
188 double r6 = r5*r;
189 double z2 = z*z;
190 double z3 = z2*z;
191 double z4 = z3*z;
192 double z5 = z4*z;
193 double z6 = z5*z;
194 double logr = log(r);
195
196 psi[0] = Bdata->psi_mult * (
197 (1-C[12]) * (r4/8)
198 + C[12] * (r2*logr/2)
199 + C[0] * (1)
200 + C[1] * (r2)
201 + C[2] * (r2*logr - z2)
202 + C[3] * (r4 - 4*r2*z2)
203 + C[4] * (3*r4*logr - 9*r2*z2 - 12*r2*logr*z2 + 2*z4)
204 + C[5] * (r6 - 12*r4*z2 + 8*r2*z4)
205 + C[6] * (8*z6 - 140*r2*z4 - 120*r2*logr*z4 + 180*r4*logr*z2
206 + 75*r4*z2 - 15*r6*logr)
207 + C[7] * (z)
208 + C[8] * (z*r2)
209 + C[9] * (z3 - 3*z*r2*logr)
210 + C[10] * (3*z*r4 - 4*z3*r2)
211 + C[11] * (8*z5 - 45*z*r4 - 80*z3*r2*logr + 60*z*r4*logr) );
212
213 return 0;
214}
215
227a5err B_GS_eval_psi_dpsi(real psi_dpsi[4], real r, real phi, real z,
228 B_GS_data* Bdata) {
229
230 // The psi value itself we can get from here
231 B_GS_eval_psi(&(psi_dpsi[0]), r, phi, z, Bdata);
232
233 /* Normalize the coordinates */
234 z -= Bdata->z0;
235 r /= Bdata->R0;
236 z /= Bdata->R0;
237
238 real* C = Bdata->psi_coeff;
239
240 /* Precalculated terms used in the components */
241 real r2 = r*r;
242 real r3 = r2*r;
243 real r4 = r3*r;
244 real r5 = r4*r;
245 real z2 = z*z;
246 real z3 = z2*z;
247 real z4 = z3*z;
248 real z5 = z4*z;
249 real logr = log(r);
250
251 psi_dpsi[1] = ( Bdata->psi_mult / (r * Bdata->R0) ) *
252 ( (1-C[12]) * (r3/2)
253 + C[12] * (r/2 + r*logr)
254 + C[1] * (2*r)
255 + C[2] * (2*r*logr + r)
256 + C[3] * (4*r3 - 8*r*z2)
257 + C[4] * (12*r3*logr + 3*r3 - 30*r*z2 - 24*r*logr*z2)
258 + C[5] * (6*r5 - 48*r3*z2 + 16*r*z4)
259 + C[6] * (-400*r*z4 -240*r*logr*z4 + 720*r3*logr*z2 + 480*r3*z2
260 -90*r5*logr - 15*r5)
261 + C[8] * (2*z*r)
262 + C[9] * (-6*z*r*logr - 3*z*r)
263 + C[10] * (12*z*r3 - 8*z3*r)
264 + C[11] * (-120*z*r3-160*z3*r*logr-80*z3*r+240*z*r3*logr) );
265
266 psi_dpsi[2] = 0;
267
268 psi_dpsi[3] =
269 ( Bdata->psi_mult / (r * Bdata->R0) ) *
270 ( C[2] * (-2*z)
271 + C[3] * (-8*r2*z)
272 + C[4] * (-18*r2*z - 24*r2*logr*z + 8*z3)
273 + C[5] * (-24*r4*z + 32*r2*z3)
274 + C[6] * (48*z5 - 560*r2*z3 - 480*r2*logr*z3 +360*r4*logr*z
275 + 150*r4*z)
276 + C[7] * (1)
277 + C[8] * (r2)
278 + C[9] * (3*z2 - 3*r2*logr)
279 + C[10] * (3*r4 - 12*z2*r2)
280 + C[11] * (40*z4 - 45*r4 - 240*z2*r2*logr + 60*r4*logr) );
281
282 return 0;
283}
284
296a5err B_GS_eval_rho_drho(real rho_drho[4], real r, real phi, real z,
297 B_GS_data* Bdata) {
298 real psi_dpsi[4];
299
300 B_GS_eval_psi_dpsi(psi_dpsi, r, phi, z, Bdata);
301
302 /* Check that the values seem valid */
303 real delta = Bdata->psi1 - Bdata->psi0;
304 if( (psi_dpsi[0] - Bdata->psi0) / delta < 0 ) {
305 return error_raise( ERR_INPUT_UNPHYSICAL, __LINE__, EF_B_GS );
306 }
307
308 rho_drho[0] = sqrt( (psi_dpsi[0] - Bdata->psi0) / delta );
309 rho_drho[1] = psi_dpsi[1] / (2*delta*rho_drho[0]);
310 rho_drho[2] = 0;
311 rho_drho[3] = psi_dpsi[3] / (2*delta*rho_drho[0]);
312
313 return 0;
314}
315
327a5err B_GS_eval_B(real B[3], real r, real phi,
328 real z, B_GS_data* Bdata) {
329 /* Normalize the coordinates */
330 z -= Bdata->z0;
331 r /= Bdata->R0;
332 z /= Bdata->R0;
333
334 real* C = Bdata->psi_coeff;
335
336 /* Precalculated terms used in the components */
337 real r2 = r*r;
338 real r3 = r2*r;
339 real r4 = r3*r;
340 real r5 = r4*r;
341 real z2 = z*z;
342 real z3 = z2*z;
343 real z4 = z3*z;
344 real z5 = z4*z;
345 real logr = log(r);
346
347 /* r field component */
348 B[0] = C[2] * (-2*z)
349 + C[3] * (-8*r2*z)
350 + C[4] * (-18*r2*z - 24*r2*logr*z + 8*z3)
351 + C[5] * (-24*r4*z + 32*r2*z3)
352 + C[6] * (48*z5 - 560*r2*z3 - 480*r2*logr*z3 +360*r4*logr*z
353 + 150*r4*z)
354 + C[7] * (1)
355 + C[8] * (r2)
356 + C[9] * (3*z2 - 3*r2*logr)
357 + C[10] * (3*r4 - 12*z2*r2)
358 + C[11] * (40*z4 - 45*r4 - 240*z2*r2*logr + 60*r4*logr);
359 B[0] *= -Bdata->psi_mult / (r * Bdata->R0 * Bdata->R0);
360
361 /* phi field component */
362 B[1] = Bdata->B_phi0 / r;
363
364 /* z field component */
365 B[2] = (1-C[12]) * (r3/2)
366 + C[12] * (r/2 + r*logr)
367 + C[1] * (2*r)
368 + C[2] * (2*r*logr + r)
369 + C[3] * (4*r3 - 8*r*z2)
370 + C[4] * (12*r3*logr + 3*r3 - 30*r*z2 - 24*r*logr*z2)
371 + C[5] * (6*r5 - 48*r3*z2 + 16*r*z4)
372 + C[6] * (-400*r*z4 -240*r*logr*z4 + 720*r3*logr*z2 + 480*r3*z2
373 -90*r5*logr - 15*r5)
374 + C[8] * (2*z*r)
375 + C[9] * (-6*z*r*logr - 3*z*r)
376 + C[10] * (12*z*r3 - 8*z3*r)
377 + C[11] * (-120*z*r3-160*z3*r*logr-80*z3*r+240*z*r3*logr);
378 B[2] *= Bdata->psi_mult / (r * Bdata->R0 * Bdata->R0);
379
380 /* Ripple */
381 if(Bdata->Nripple > 0) {
382 r *= Bdata->R0;
383 z *= Bdata->R0;
384 real radius = sqrt( ( r - Bdata->R0 ) * ( r - Bdata->R0 )
385 + ( z - Bdata->z0 ) * ( z - Bdata->z0 ));
386 real theta = atan2( z - Bdata->z0, r - Bdata->R0 );
387 real delta = Bdata->delta0 * exp(-0.5*theta*theta)
388 * pow( radius / Bdata->a0, Bdata->alpha0 );
389 B[1] = B[1] * ( 1 + delta * cos(Bdata->Nripple * phi) );
390 }
391
392 return 0;
393}
394
406a5err B_GS_eval_B_dB(real B_dB[12], real r, real phi, real z,
407 B_GS_data* Bdata) {
408
409 real* C = Bdata->psi_coeff;
410
411 z -= Bdata->z0;
412 real R0 = Bdata->R0;
413 real z0 = Bdata->z0;
414 real B_phi0 = Bdata->B_phi0;
415 real psi_mult = Bdata->psi_mult;
416
417 /* Precalculated terms used in the components */
418 real r2 = r*r;
419 real r3 = r2*r;
420 real r4 = r3*r;
421 real r5 = r4*r;
422 real z2 = z*z;
423 real z3 = z2*z;
424 real z4 = z3*z;
425 real z5 = z4*z;
426 real logr = log(r/R0);
427 real R02 = R0*R0;
428 real R03 = R02*R0;
429 real R04 = R03*R0;
430 real R05 = R04*R0;
431 real R06 = R05*R0;
432
433 /* r field component */
434 real B0 = C[2] * (-2*z) / R02
435 + C[3] * (-8*r2*z) / R04
436 + C[4] * (-18*r2*z - 24*r2*logr*z + 8*z3) / R04
437 + C[5] * (-24*r4*z + 32*r2*z3) / R06
438 + C[6] * (48*z5 - 560*r2*z3 - 480*r2*logr*z3 +360*r4*logr*z
439 + 150*r4*z) / R06
440 + C[7] * (1) / R0
441 + C[8] * (r2) / R03
442 + C[9] * (3*z2 - 3*r2*logr) / R03
443 + C[10] * (3*r4 - 12*z2*r2) / R05
444 + C[11] * (40*z4 - 45*r4 - 240*z2*r2*logr + 60*r4*logr) / R05;
445 B0 *= -psi_mult / r;
446 real B1 = C[3] * (-16*r*z) / R04
447 + C[4] * (-60*r*z - 48*r*logr*z) / R04
448 + C[5] * (-96*r3*z + 64*r*z3) / R06
449 + C[6] * (-1600*r*z3 - 960*r*logr*z3 +1440*r3*logr*z
450 + 960*r3*z) / R06
451 + C[8] * (2*r) / R03
452 + C[9] * (-6*r*logr - 3*r) / R03
453 + C[10] * (12*r3 - 24*z2*r) / R05
454 + C[11] * (-120*r3 - 480*z2*r*logr -240*z2*r +240*r3*logr)/R05;
455 B1 = -B0 / r - B1 * psi_mult / r;
456 real B2 = 0;
457 real B3 = C[2] * (-2) / R02
458 + C[3] * (-8*r2) / R04
459 + C[4] * (-18*r2 - 24*r2*logr + 24*z2) / R04
460 + C[5] * (-24*r4 + 96*r2*z2) / R06
461 + C[6] * (240*z4 - 1680*r2*z2 - 1440*r2*logr*z2 + 360*r4*logr
462 + 150*r4) / R06
463 + C[9] * (6*z) / R03
464 + C[10] * (-24*z*r2) / R05
465 + C[11] * (160*z3 - 480*z*r2*logr) / R05;
466 B3 *= -psi_mult / r;
467
468 /* phi field component */
469 real B4 = B_phi0 * R0 / r;
470 real B5 = -B_phi0 * R0 / r2;
471 real B6 = 0;
472 real B7 = 0;
473
474 /* z field component */
475 real B8 = (1-C[12]) * (r3/2) / R04
476 + C[12] * (r/2 + r*logr) / R02
477 + C[1] * (2*r) / R02
478 + C[2] * (2*r*logr + r) / R02
479 + C[3] * (4*r3 - 8*r*z2) / R04
480 + C[4] * (12*r3*logr + 3*r3 - 30*r*z2 - 24*r*logr*z2) / R04
481 + C[5] * (6*r5 - 48*r3*z2 + 16*r*z4) / R06
482 + C[6] * (-400*r*z4 -240*r*logr*z4 + 720*r3*logr*z2 + 480*r3*z2
483 -90*r5*logr - 15*r5) / R06
484 + C[8] * (2*z*r) / R03
485 + C[9] * (-6*z*r*logr - 3*z*r) / R03
486 + C[10] * (12*z*r3 - 8*z3*r) / R05
487 + C[11] * (-120*z*r3-160*z3*r*logr-80*z3*r+240*z*r3*logr) / R05;
488 B8 *= psi_mult / r;
489 real B9 = (1-C[12]) * (3*r2/2) / R04
490 + C[12] * (1.5 + logr) / R02
491 + C[1] * (2) / R02
492 + C[2] * (2*logr + 3) / R02
493 + C[3] * (12*r2 - 8*z2) / R04
494 + C[4] * (36*r2*logr + 21*r2 - 54*z2 - 24*logr*z2) / R04
495 + C[5] * (30*r4 - 144*r2*z2 + 16*z4) / R06
496 + C[6] * (-640*z4 - 240*logr*z4 + 2160*r2*logr*z2 + 2160*r2*z2
497 -450*r4*logr -165*r4) / R06
498 + C[8] * (2*z) / R03
499 + C[9] * (-6*z*logr - 9*z) / R03
500 + C[10] * (36*z*r2 - 8*z3) / R05
501 + C[11] * (-120*z*r2 - 160*z3*logr -240*z3 +720*z*r2*logr)/R05;
502 B9 = B9 * psi_mult / r - B8 / r;
503 real B10 = 0;
504 real B11 = -B1 - B0 / r;
505
506 /* Ripple */
507 if(Bdata->Nripple > 0) {
508 real radius = sqrt( ( r - R0 ) * ( r - R0 )
509 + ( z - z0 ) * ( z - z0 ));
510 real theta = atan2( z - z0, r - R0 );
511 real delta = Bdata->delta0 * exp(-0.5*theta*theta)
512 * pow( radius / Bdata->a0, Bdata->alpha0 );
513
514 real Bphi = B4;
515 real Bpert = Bphi * delta * cos(Bdata->Nripple * phi);
516 B4 += Bpert;
517 B6 += - Bphi * delta * Bdata->Nripple * sin(Bdata->Nripple * phi);
518
519 real dBpertdR = Bpert * (
520 ( (r - R0) /radius) * ( Bdata->alpha0 / radius )
521 + ( (z - z0) /(radius*radius) ) * theta
522 );
523
524 real dBpertdz = Bpert * (
525 ( (z - z0) /radius) * ( Bdata->alpha0 / radius )
526 - ( (r - R0) /(radius*radius) ) * theta
527 );
528
529 B5 += B5 * Bpert / Bphi + dBpertdR;
530 B7 += dBpertdz;
531 }
532
533 B_dB[0] = B0;
534 B_dB[1] = B1;
535 B_dB[2] = B2;
536 B_dB[3] = B3;
537 B_dB[4] = B4;
538 B_dB[5] = B5;
539 B_dB[6] = B6;
540 B_dB[7] = B7;
541 B_dB[8] = B8;
542 B_dB[9] = B9;
543 B_dB[10] = B10;
544 B_dB[11] = B11;
545
546 return 0;
547}
548
557a5err B_GS_get_axis_rz(real rz[2], B_GS_data* Bdata) {
558 a5err err = 0;
559 rz[0] = Bdata->raxis;
560 rz[1] = Bdata->zaxis;
561 return err;
562}
B_GS_init
void B_GS_init(B_GS_data *Bdata, B_GS_offload_data *offload_data, real *offload_array)
Initialize magnetic field data struct on target.
Definition B_GS.c:144
B_GS_eval_psi
a5err B_GS_eval_psi(real *psi, real r, real phi, real z, B_GS_data *Bdata)
Evaluate poloidal flux psi.
Definition B_GS.c:174
B_GS_eval_B
a5err B_GS_eval_B(real B[3], real r, real phi, real z, B_GS_data *Bdata)
Evaluate magnetic field.
Definition B_GS.c:327
B_GS_eval_B_dB
a5err B_GS_eval_B_dB(real B_dB[12], real r, real phi, real z, B_GS_data *Bdata)
Evaluate magnetic field and its derivatives.
Definition B_GS.c:406
B_GS_get_axis_rz
a5err B_GS_get_axis_rz(real rz[2], B_GS_data *Bdata)
Return magnetic axis R-coordinate.
Definition B_GS.c:557
B_GS_free_offload
void B_GS_free_offload(B_GS_offload_data *offload_data, real **offload_array)
Free offload array.
Definition B_GS.c:131
B_GS_eval_rho_drho
a5err B_GS_eval_rho_drho(real rho_drho[4], real r, real phi, real z, B_GS_data *Bdata)
Evaluate normalized poloidal flux rho and its derivatives.
Definition B_GS.c:296
B_GS_eval_psi_dpsi
a5err B_GS_eval_psi_dpsi(real psi_dpsi[4], real r, real phi, real z, B_GS_data *Bdata)
Evaluate poloidal flux psi and its derivatives.
Definition B_GS.c:227
B_GS_init_offload
int B_GS_init_offload(B_GS_offload_data *offload_data, real **offload_array)
Initialize magnetic field offload data.
Definition B_GS.c:91
B_GS.h
Header file for B_GS.c.
ascot5.h
Main header file for ASCOT5.
real
double real
Definition ascot5.h:85
consts.h
Header file containing physical and mathematical constants.
error.h
Error module for ASCOT5.
a5err
unsigned long int a5err
Simulation error flag.
Definition error.h:17
EF_B_GS
@ EF_B_GS
Definition error.h:39
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
math.h
Header file for math.c.
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_GS_data
Analytic magnetic field parameters on the target.
Definition B_GS.h:35
B_GS_offload_data
Analytic magnetic field parameters that will be offloaded to target.
Definition B_GS.h:14
B_GS_offload_data::R0
real R0
Definition B_GS.h:15
B_GS_offload_data::zaxis
real zaxis
Definition B_GS.h:18
B_GS_offload_data::raxis
real raxis
Definition B_GS.h:17
B_GS_offload_data::a0
real a0
Definition B_GS.h:26
B_GS_offload_data::Nripple
int Nripple
Definition B_GS.h:25
B_GS_offload_data::psi0
real psi0
Definition B_GS.h:20
B_GS_offload_data::B_phi0
real B_phi0
Definition B_GS.h:19
B_GS_offload_data::delta0
real delta0
Definition B_GS.h:28
B_GS_offload_data::psi_coeff
real psi_coeff[13]
Definition B_GS.h:23
B_GS_offload_data::psi1
real psi1
Definition B_GS.h:21
B_GS_offload_data::psi_mult
real psi_mult
Definition B_GS.h:22
B_GS_offload_data::offload_array_length
int offload_array_length
Definition B_GS.h:29
B_GS_offload_data::z0
real z0
Definition B_GS.h:16
B_GS_offload_data::alpha0
real alpha0
Definition B_GS.h:27
Generated on Wed Oct 9 2024 17:13:56 for ASCOT5 by doxygen 1.10.0