forked from mabarnes/stella
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathsfincs_interface.fpp
1580 lines (1305 loc) · 55.6 KB
/
sfincs_interface.fpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
module sfincs_interface
implicit none
public :: get_neo_from_sfincs
private
# ifdef USE_SFINCS
integer :: nproc_sfincs
integer :: irad_min, irad_max
real :: Er_window
logical :: includeXDotTerm
logical :: includeElectricFieldTermInXiDot
! logical :: includeRadialExBDrive
integer :: magneticDriftScheme
logical :: includePhi1
logical :: includePhi1InKineticEquation
! logical :: includePhi1InCollisionOperator
integer :: geometryScheme
integer :: VMECRadialOption
integer :: coordinateSystem
integer :: inputRadialCoordinate
integer :: inputRadialCoordinateForGradients
character (200) :: equilibriumFile
real :: aHat, psiAHat, Delta
real :: nu_n, dPhiHatdrN
integer :: nxi, nx, ntheta, nzeta
logical :: calculate_radial_electric_field
logical :: read_sfincs_output_from_file
logical :: sfincs_finished = .true.
real, dimension (:), allocatable :: fprim_local, tprim_local
# endif
contains
subroutine get_neo_from_sfincs (nradii, drho, f_neoclassical, phi_neoclassical, &
dfneo_dalpha, dphineo_dalpha)
# ifdef USE_SFINCS
use mp, only: proc0, iproc
use mp, only: comm_split, comm_free
use stella_geometry, only: geo_surf
use species, only: spec, nspec
# endif
use mp, only: mp_abort
use zgrid, only: nzgrid
implicit none
integer, intent (in) :: nradii
real, intent (in) :: drho
real, dimension (:,-nzgrid:,:,:,:,-nradii/2:), intent (out) :: f_neoclassical
real, dimension (:,-nzgrid:,-nradii/2:), intent (out) :: phi_neoclassical
real, dimension (:,-nzgrid:,:,:,:), intent (out) :: dfneo_dalpha
real, dimension (:,-nzgrid:), intent (out) :: dphineo_dalpha
# ifdef USE_SFINCS
integer :: sfincs_comm
integer :: color, ierr
integer :: irad
real :: dPhiHatdrN_best_guess
logical :: Er_converged
integer :: nsfincs_calls
if (.not.allocated(fprim_local)) allocate (fprim_local(nspec))
if (.not.allocated(tprim_local)) allocate (tprim_local(nspec))
if (proc0) call read_sfincs_parameters (nradii)
call broadcast_sfincs_parameters
if (iproc < nproc_sfincs) then
color = 0
else
color = 1
end if
call comm_split (color, sfincs_comm, ierr)
if (iproc < nproc_sfincs) then
do irad = irad_min, irad_max
! get local values of -dlog(ns)/drho and -dlog(Ts)/drho
! using dlog(n)/drho = dlog(n0)/drho + delrho*d/drho(dlog(n)/drho)
fprim_local = 1.0/geo_surf%drhotordrho*(spec%fprim &
+ irad*drho*(spec%fprim**2-spec%d2ndr2)/geo_surf%drhotordrho)
tprim_local = 1.0/geo_surf%drhotordrho*(spec%tprim &
+ irad*drho*(spec%tprim**2-spec%d2Tdr2)/geo_surf%drhotordrho)
! fprim_local = 1.0/geo_surf%drhotordrho*(spec%fprim - irad*drho*spec%d2ndr2)
! tprim_local = 1.0/geo_surf%drhotordrho*(spec%tprim - irad*drho*spec%d2Tdr2)
if (calculate_radial_electric_field) then
! get best guess at radial electric field
! using force balance with radial pressure gradient
if (dPhiHatdrN > -9999.0) then
dPhiHatdrN_best_guess = dPhiHatdrN
else
dPhiHatdrN_best_guess = fprim_local(1)+tprim_local(1)
end if
call iterate_sfincs_until_electric_field_converged (sfincs_comm, &
irad, drho, irad_max, dPhiHatdrN_best_guess, &
Er_converged, nsfincs_calls)
if (proc0) then
write (*,*)
write (*,*) 'At irad= ', irad, 'Er_converged= ', Er_converged, &
'nsfincs_calls_required= ', nsfincs_calls, 'dPhiHatdrN= ', dPhiHatdrN
write (*,*)
end if
! write_and_finish_sfincs manipulates sfincs output
! to get the neoclassical distribution function and potential
! on the stella (zed,alpha,vpa,mu) grid; it then
! deallocates sfincs arrays, etc. to make it ready
! for running again if need be
call write_and_finish_sfincs (f_neoclassical(:,:,:,:,:,irad), &
phi_neoclassical(:,:,irad), dfneo_dalpha, dphineo_dalpha, irad)
else
! init_and_run_sfincs initializes sfincs,
! including passing geometry info if necessary;
! and runs sfincs (if requested)
call init_and_run_sfincs (sfincs_comm, irad, drho, irad_max)
! write_and_finish_sfincs manipulates sfincs output
! to get the neoclassical distribution function and potential
! on the stella (zed,alpha,vpa,mu) grid; it then
! deallocates sfincs arrays, etc. to make it ready
! for running again if need be
call write_and_finish_sfincs (f_neoclassical(:,:,:,:,:,irad), &
phi_neoclassical(:,:,irad), dfneo_dalpha, dphineo_dalpha, irad)
end if
end do
end if
call comm_free (sfincs_comm, ierr)
! NB: NEED TO CHECK THIS BROADCAST OF SFINCS RESULTS
do irad = irad_min, irad_max
call broadcast_sfincs_output &
(f_neoclassical(:,:,:,:,:,irad), phi_neoclassical(:,:,irad))
end do
call broadcast_sfincs_output (dfneo_dalpha, dphineo_dalpha)
deallocate (fprim_local, tprim_local)
if (proc0) then
write (*,*)
write (*,*) 'maxval(fneo): ', maxval(f_neoclassical(:,:,:,:,:,irad_min:irad_max)), &
'maxval(phineo): ', maxval(phi_neoclassical(:,:,irad_min:irad_max))
write (*,*)
end if
if ((irad_min /= -nradii/2) .or. (irad_max /= nradii/2)) &
call mp_abort ('WARNING: irad_min must equal -nradii/2 and irad_max must equal &
& nradii/2 to proceed to stella calculation. aborting.')
# else
real :: dum
! this pointless dum assignment only here to avoid
! annoying warning messages during compilation
! about unused variable drho
dum = drho
f_neoclassical = 0. ; phi_neoclassical = 0.
dfneo_dalpha = 0. ; dphineo_dalpha = 0.
call mp_abort ("to run with include_neoclassical_terms=.true., &
& USE_SFINCS must be defined at compilation time. Aborting.")
# endif
end subroutine get_neo_from_sfincs
# ifdef USE_SFINCS
subroutine iterate_sfincs_until_electric_field_converged (sfincs_comm, irad, drho, &
nrad_max, dPhiHatdrN_best_guess, dphiHatdrN_is_converged, number_of_sfincs_calls_for_convergence)
use mp, only: proc0, iproc
implicit none
integer, intent (in) :: sfincs_comm, irad, nrad_max
real, intent (in) :: drho, dPhiHatdrN_best_guess
logical, intent (out) :: dPhiHatdrN_is_converged
integer, intent (out) :: number_of_sfincs_calls_for_convergence
integer :: itmax_bracket = 10
integer :: itmax_root = 10
real :: window = 0.3
real :: tol = 0.1
integer :: it
real :: a,b,c,d,e,fa,fb,fc,p,q,r,s,tol1,xm,eps
real :: converged_dPhiHatdrN
dPhiHatdrN_is_converged = .false.
a = dPhiHatdrN_best_guess*(1.0-Er_window)
b = dPhiHatdrN_best_guess*(1.0+Er_window)
! initialize sfincs, run it, and return the total charge flux as fa
call get_total_charge_flux (sfincs_comm, irad, drho, nrad_max, a, fa)
call get_total_charge_flux (sfincs_comm, irad, drho, nrad_max, b, fb)
number_of_sfincs_calls_for_convergence = 2
do it = 1, itmax_bracket
eps = epsilon(a)
if ((fa > 0.0 .and. fb > 0.0) .or. (fa < 0.0 .and. fb < 0.0)) then
if (proc0) then
write (*,*)
write (*,*) 'dPhiHatdrN values ', a, ' and ', b, ' do not bracket root.'
write (*,*) 'flux at ', a, ' is ', fa, '.'
write (*,*) 'flux at ', b, ' is ', fb, '.'
end if
a = a*(1.0-Er_window)
b = b*(1.0+Er_window)
if (proc0) then
write (*,*) 'Trying again with values ', a, ' and ', b, ' .'
write (*,*)
end if
call get_total_charge_flux (sfincs_comm, irad, drho, nrad_max, a, fa)
call get_total_charge_flux (sfincs_comm, irad, drho, nrad_max, b, fb)
number_of_sfincs_calls_for_convergence = number_of_sfincs_calls_for_convergence + 2
! ! eliminate the endpoint corresonding to the flux that is furthest from zero in magnitude
! if (abs(fa) > abs(fb)) then
! ! keep b as an endpoint and eliminate a
! a = b ; fa = fb
! b = a*(1.0+window)
! if (proc0) then
! write (*,*) 'Trying again with values ', a, ' and ', b, ' .'
! write (*,*)
! end if
! call get_total_charge_flux (sfincs_comm, irad, drho, nrad_max, b, fb)
! else
! ! keep a as an endpoint and eliminate b
! b = a ; fb = fa
! a = b*(1.0-window)
! if (proc0) then
! write (*,*) 'Trying again with values ', a, ' and ', b, ' .'
! write (*,*)
! end if
! call get_total_charge_flux (sfincs_comm, irad, drho, nrad_max, a, fa)
! end if
else
exit
end if
end do
c=b
fc=fb
do it = 1, itmax_root
if ((fb > 0.0 .and. fc > 0.0) .or. (fb < 0.0 .and. fc < 0.0)) then
c=a
fc=fa
d=b-a
e=d
end if
if (abs(fc) < abs(fb)) then
a=b
b=c
c=a
fa=fb
fb=fc
fc=fa
end if
tol1=2.0*eps*abs(b)+0.5*tol
xm=0.5*(c-b)
if (abs(xm) <= tol1 .or. fb == 0.0) then
converged_dPhiHatdrN = b
dPhiHatdrN_is_converged = .true.
! number_of_sfincs_calls_for_convergence = it+1
exit
end if
if (abs(e) >= tol1 .and. abs(fa) > abs(fb)) then
s=fb/fa
if (a==c) then
p=2.0*xm*s
q=1.0-s
else
q=fa/fc
r=fb/fc
p=s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0))
q=(q-1.0)*(r-1.0)*(s-1.0)
end if
if (p > 0.0) q=-q
p=abs(p)
if (2.0*p < min(3.0*xm*q-abs(tol1*q),abs(e*q))) then
e=d
d=p/q
else
d=xm
e=d
end if
else
d=xm
e=d
end if
a=b
fa=fb
b=b+merge(d,sign(tol1,xm), abs(d) > tol1)
call get_total_charge_flux (sfincs_comm, irad, drho, nrad_max, b, fb)
number_of_sfincs_calls_for_convergence = number_of_sfincs_calls_for_convergence + 1
end do
end subroutine iterate_sfincs_until_electric_field_converged
subroutine get_total_charge_flux (sfincs_comm, irad, drho, nrad_max, &
dPhiHatdrN_in, total_charge_flux)
use mp, only: iproc, broadcast_with_comm
use sfincs_main, only: finish_sfincs
use globalVariables, only: Zs, particleFlux_vd_psiHat
use species, only: nspec
implicit none
integer, intent (in) :: sfincs_comm, irad, nrad_max
real, intent (in) :: drho
real, intent (in) :: dPhiHatdrN_in
real, intent (out) :: total_charge_flux
dPhiHatdrN = dPhiHatdrN_in
call init_and_run_sfincs (sfincs_comm, irad, drho, nrad_max)
call broadcast_with_comm (Zs, sfincs_comm)
call broadcast_with_comm (particleFlux_vd_psiHat, sfincs_comm)
total_charge_flux = sum(Zs(:nspec)*particleFlux_vd_psiHat(:nspec))
end subroutine get_total_charge_flux
subroutine init_and_run_sfincs (sfincs_comm, irad, drho, nrad_max)
use mp, only: proc0, iproc
use sfincs_main, only: init_sfincs, prepare_sfincs, run_sfincs, finish_sfincs
use globalVariables, only: Zs, particleFlux_vd_psiHat
use species, only: nspec
implicit none
integer, intent (in) :: sfincs_comm, irad, nrad_max
real, intent (in) :: drho
if (.not.sfincs_finished) then
call finish_sfincs
sfincs_finished = .true.
end if
call init_sfincs (sfincs_comm)
call pass_inputoptions_to_sfincs (irad*drho)
call pass_outputoptions_to_sfincs
call prepare_sfincs
! if geometryScheme = 5, then sfincs will read in equilibrium
! parameters from vmec file separately
! otherwise, assume system is axisymmetric and pass geometry
! from stella (miller local equilibrium or similar)
if (geometryScheme /= 5) call pass_geometry_to_sfincs (irad*drho)
if (read_sfincs_output_from_file) then
if (proc0) call read_sfincs_output (irad, nrad_max)
else
if (proc0) then
write (*,*)
write (*,*) 'Running sfincs at irad= ', irad, ', with dPhiHatdrN= ', dPhiHatdrN
end if
call run_sfincs
if (proc0) then
write (*,*) 'sfincs finished running. total charge flux= ', sum(Zs(:nspec)*particleFlux_vd_psiHat(:nspec))
write (*,*)
end if
! write Phi1Hat and delta_f to file
! so we have the option of using it
! again without re-running sfincs
if (proc0) call write_sfincs (irad, nrad_max)
end if
sfincs_finished = .false.
end subroutine init_and_run_sfincs
subroutine write_and_finish_sfincs (fneo, phineo, dfneo, dphineo, irad)
use mp, only: proc0
use sfincs_main, only: finish_sfincs
use zgrid, only: nzgrid
implicit none
real, dimension (:,-nzgrid:,:,:,:), intent (out) :: fneo
real, dimension (:,-nzgrid:), intent (out) :: phineo
real, dimension (:,-nzgrid:,:,:,:), intent (in out) :: dfneo
real, dimension (:,-nzgrid:), intent (in out) :: dphineo
integer, intent (in) :: irad
if (proc0) then
! only need to compute dfneo_dalpha and dphineo_dalpha
! for central radius and for stellarator calculation
if (irad == 0 .and. geometryScheme == 5) then
call get_sfincs_output (fneo, phineo, dfneo, dphineo)
else
call get_sfincs_output (fneo, phineo)
end if
end if
call finish_sfincs
end subroutine write_and_finish_sfincs
subroutine read_sfincs_parameters (nradii)
use constants, only: pi
use mp, only: nproc
use file_utils, only: input_unit_exist
use species, only: nspec
use physics_parameters, only: rhostar, vnew_ref
use stella_geometry, only: geo_surf, aref, bref
implicit none
integer, intent (in) :: nradii
namelist /sfincs_input/ nproc_sfincs, &
calculate_radial_electric_field, &
includeXDotTerm, &
includeElectricFieldTermInXiDot, &
irad_min, irad_max, &
magneticDriftScheme, &
includePhi1, &
includePhi1InKineticEquation, &
! includePhi1InCollisionOperator, &
geometryScheme, &
VMECRadialOption, &
equilibriumFile, &
coordinateSystem, &
inputRadialCoordinate, &
inputRadialCoordinateForGradients, &
aHat, psiAHat, nu_N, nxi, nx, Delta, &
dPhiHatdrN, &
ntheta, nzeta, &
read_sfincs_output_from_file, Er_window
logical :: exist
integer :: in_file
! if read_sfincs_output_from_file=.true.,
! will try to read in Phi1Hat and delta_f
! from pre-saved file named sfincs.output
! otherwise, run sfincs to compute these
! quantities on the fly
read_sfincs_output_from_file = .false.
! number of processors to use for sfincs calculation
nproc_sfincs = 1
! minimum and maximum radial index (irad=0 corresponds to central radius)
irad_min = -nradii/2 ; irad_max = nradii/2
! if calculate_radial_electric_field, then
! will scan in radial electric field to find value
! for which ambipolarity is satisfied, and then
! use this value to obtain neoclassical fluxes,
! distribution function, and potential
calculate_radial_electric_field = .true.
! do not include radial electric field term if set to .false.
includeXDotTerm = .true.
includeElectricFieldTermInXiDot = .true.
! include v_E . grad r term
! includeRadialExBDrive = .true.
! no poloidal or toroidal magnetic drifts
magneticDriftScheme = 0
! combo of next two variables means
! phi1 will be calculated via quasineutrality
includePhi1 = .true.
includePhi1InKineticEquation = .false.
! includePhi1InCollisionOperator = .false.
! will be overridden by direct input of geometric quantities
! unless geometryScheme = 5 (vmec equilibrium)
geometryScheme = 1
! only relevant if geometryScheme = 5
! radial option to use for vmec equilibrium
! 0 corresponds to using radial interpolation to get desired surface
! 1 corresponds to using nearest surface on VMEC HALF grid
! 2 corresponds to using nearest surface on VMEC FULL grid
! should not change this unless self-consistently change in the
! vmec input namelist
VMECRadialOption = 0
! path of vmec equilibrium file
equilibriumFile = 'wout_161s1.nc'
! seems to be a nonsensical option
coordinateSystem = 3
! option 3 corresponds to using sqrt of toroidal flux
! normalized by toroidal flux enclosed by the LCFS
inputRadialCoordinate = 3
! option 3 corresponds to same choice
! when calculating gradients of density, temperature, and potential
inputRadialCoordinateForGradients = 3
! corresponds to r_LCFS as reference length in sfincs
! only used in sfincs when geometryScheme=1
aHat = 1.0
! psitor_LCFS / (B_ref * a_ref^2)
psiAHat = geo_surf%psitor_lcfs
! Delta is rho* = mref*vt_ref/(e*Bref*aref), with reference
! quantities given in SI units
! unless geometryScheme = 5, in which case Bref=1T
! and aref = 1m (these are hardwired in sfincs)
! set negative to allow check later to see if any value given in input file
Delta = -1.0
! nu_n = nu_ref * aref/vt_ref
! nu_ref = 4*sqrt(2*pi)*nref*e**4*loglam/(3*sqrt(mref)*Tref**3/2)
! (with nref, Tref, and mref in Gaussian units)
! set negative to allow check later to see if any value given in input file
nu_n = -1.0
! radial derivative of normalized phi
dPhiHatdrN = -9999.9
Er_window = 0.3
! number of spectral coefficients in pitch angle
nxi = 48
! number of speeds
nx = 12
! number of poloidal angles
Ntheta = 65
! number of toroidal angles, 1 is appropriate for tokamak
Nzeta = 1
in_file = input_unit_exist("sfincs_input", exist)
if (exist) read (unit=in_file, nml=sfincs_input)
if (nproc_sfincs > nproc) then
write (*,*) 'requested number of processors for sfincs is greater &
& than total processor count.'
write (*,*) 'allocating ', nproc, ' processors for sfincs.'
end if
if (Delta < 0.0) then
Delta = rhostar
! if geometryScheme=5, Bref=1T and aref=1m are hard-wired in sfincs
! but these are not the values used in stella to define rhostar
if (geometryScheme == 5) Delta = rhostar*bref*aref
end if
if (nu_n < 0.0) then
nu_n = vnew_ref*(4./(3.*sqrt(pi)))
! if geometryScheme=5, aref=1m is hard-wired in sfincs
! but this is not the value used in stella
if (geometryScheme == 5) nu_n = nu_n/aref
end if
! FLAG -- NOT YET SURE IF THIS SHOULD BE HERE
! if (nspec == 1 .and. includePhi1) then
! write (*,*) 'includePhi1 = .true. is incompatible with a single-species run.'
! write (*,*) 'forcing includePhi1 = .false.'
! includePhi1 = .false.
! end if
! ensure that ntheta and nzeta are odd for SFINCS
ntheta = 2*(ntheta/2)+1
nzeta = 2*(nzeta/2)+1
end subroutine read_sfincs_parameters
subroutine broadcast_sfincs_parameters
use mp, only: broadcast
use physics_parameters, only: rhostar
implicit none
call broadcast (read_sfincs_output_from_file)
call broadcast (nproc_sfincs)
call broadcast (irad_min)
call broadcast (irad_max)
call broadcast (calculate_radial_electric_field)
call broadcast (includeXDotTerm)
call broadcast (includeElectricFieldTermInXiDot)
! call broadcast (includeRadialExBDrive)
call broadcast (magneticDriftScheme)
call broadcast (includePhi1)
call broadcast (includePhi1InKineticEquation)
! call broadcast (includePhi1InCollisionOperator)
call broadcast (geometryScheme)
call broadcast (VMECRadialOption)
call broadcast (equilibriumFile)
call broadcast (coordinateSystem)
call broadcast (inputRadialCoordinate)
call broadcast (inputRadialCoordinateForGradients)
call broadcast (aHat)
call broadcast (psiAHat)
call broadcast (Delta)
call broadcast (nu_N)
call broadcast (dPhiHatdrN)
call broadcast (Er_window)
call broadcast (nxi)
call broadcast (nx)
call broadcast (ntheta)
call broadcast (nzeta)
write (*,*) 'Delta stella', Delta, rhostar
end subroutine broadcast_sfincs_parameters
subroutine pass_inputoptions_to_sfincs (delrho)
use mp, only: mp_abort
use stella_geometry, only: geo_surf
use species, only: spec, nspec
use zgrid, only: nzed
use physics_parameters, only: nine, tite
use globalVariables, only: includeXDotTerm_sfincs => includeXDotTerm
use globalVariables, only: includeElectricFieldTermInXiDot_sfincs => includeElectricFieldTermInXiDot
! use globalVariables, only: includeRadialExBDrive_sfincs => includeRadialExBDrive
use globalVariables, only: magneticDriftScheme_sfincs => magneticDriftScheme
use globalVariables, only: includePhi1_sfincs => includePhi1
use globalVariables, only: includePhi1InKineticEquation_sfincs => includePhi1InKineticEquation
! use globalVariables, only: includePhi1InCollisionOperator_sfincs => includePhi1InCollisionOperator
use globalVariables, only: geometryScheme_sfincs => geometryScheme
use globalVariables, only: equilibriumFile_sfincs => equilibriumFile
use globalVariables, only: VMECRadialOption_sfincs => VMECRadialOption
use globalVariables, only: coordinateSystem_sfincs => coordinateSystem
use globalVariables, only: RadialCoordinate => inputRadialCoordinate
use globalVariables, only: RadialCoordinateForGradients => inputRadialCoordinateForGradients
use globalVariables, only: rN_wish
use globalVariables, only: Nspecies, nHats, THats, MHats, Zs
use globalVariables, only: adiabaticNHat, adiabaticTHat, adiabaticZ
use globalVariables, only: nxi_sfincs => Nxi
use globalVariables, only: nx_sfincs => Nx
use globalVariables, only: ntheta_sfincs => Ntheta
use globalVariables, only: nzeta_sfincs => Nzeta
use globalVariables, only: dnHatdrNs, dTHatdrNs
use globalVariables, only: aHat_sfincs => aHat
use globalVariables, only: psiAHat_sfincs => psiAHat
use globalVariables, only: Delta_sfincs => Delta
use globalVariables, only: nu_n_sfincs => nu_n
use globalVariables, only: dPhiHatdrN_sfincs => dPhiHatdrN
use globalVariables, only: withAdiabatic
implicit none
real, intent (in) :: delrho
includeXDotTerm_sfincs = includeXDotTerm
includeElectricFieldTermInXiDot_sfincs = includeElectricFieldTermInXiDot
! includeRadialExBDrive_sfincs = includeRadialExBDrive
magneticDriftScheme_sfincs = magneticDriftScheme
includePhi1_sfincs = includePhi1
includePhi1InKineticEquation_sfincs = includePhi1InKineticEquation
! includePhi1InCollisionOperator_sfincs = includePhi1InCollisionOperator
geometryScheme_sfincs = geometryScheme
VMECRadialOption_sfincs = VMECRadialOption
equilibriumFile_sfincs = trim(equilibriumFile)
coordinateSystem_sfincs = coordinateSystem
RadialCoordinate = inputRadialCoordinate
RadialCoordinateForGradients = inputRadialCoordinateForGradients
Nspecies = nspec
nHats(:nspec) = spec%dens*(1.0-delrho*spec%fprim)
THats(:nspec) = spec%temp*(1.0-delrho*spec%tprim)
mHats(:nspec) = spec%mass
Zs(:nspec) = spec%z
nzeta_sfincs = nzeta
ntheta_sfincs = ntheta
nx_sfincs = nx
nxi_sfincs = nxi
aHat_sfincs = aHat
psiAHat_sfincs = psiAHat
Delta_sfincs = Delta
nu_n_sfincs = nu_n
dPhiHatdrN_sfincs = dPhiHatdrN
if (nspec == 1) then
withAdiabatic = .true.
adiabaticNHat = nHats(1)/nine
adiabaticTHat = THats(1)/tite
adiabaticZ = -1
end if
if (inputRadialCoordinate == 3) then
rN_wish = geo_surf%rhotor + delrho*geo_surf%drhotordrho
else
call mp_abort ('only inputRadialCoordinate=3 currently supported. aborting.')
end if
if (inputRadialCoordinateForGradients == 3) then
! radial density gradient with respect to rhotor = sqrt(psitor/psitor_LCFS)
! normalized by reference density (not species density)
dnHatdrNs(:nspec) = -spec%dens*fprim_local
! radial temperature gradient with respect to rhotor = sqrt(psitor/psitor_LCFS)
! normalized by reference tmperatures (not species temperature)
dTHatdrNs(:nspec) = -spec%temp*tprim_local
else
call mp_abort ('only inputRadialCoordinateForGradients=3 currently supported. aborting.')
end if
end subroutine pass_inputoptions_to_sfincs
subroutine pass_outputoptions_to_sfincs
use export_f, only: export_f_theta_option
use export_f, only: export_f_zeta_option
use export_f, only: export_f_xi_option
use export_f, only: export_f_x_option
use export_f, only: export_delta_f, export_full_f
use export_f, only: export_f_stella
implicit none
export_f_theta_option = 0
export_f_zeta_option = 0
export_f_xi_option = 0
export_f_x_option = 0
export_delta_f = .true.
export_full_f = .false.
export_f_stella = .true.
end subroutine pass_outputoptions_to_sfincs
! if this subroutine is being called, then
! using sfincs in tokamak geometry
! so zed in stella is theta
subroutine pass_geometry_to_sfincs (delrho)
use constants, only: pi
use splines, only: linear_interp_periodic
use zgrid, only: nz2pi, zed
use stella_geometry, only: bmag, dbdzed, gradpar
use stella_geometry, only: dBdrho, d2Bdrdth, dgradpardrho, dIdrho
use stella_geometry, only: geo_surf
use globalVariables, only: BHat
use globalVariables, only: dBHatdtheta
use globalVariables, only: iota
use globalVariables, only: DHat
use globalVariables, only: BHat_sup_theta
use globalVariables, only: BHat_sub_zeta
use export_f, only: export_f_theta
implicit none
real, intent (in) :: delrho
integer :: nzeta = 1
integer :: nzpi
real :: q_local
real, dimension (:), allocatable :: B_local, dBdz_local, gradpar_local
real, dimension (:), allocatable :: zed_stella
real, dimension (:), allocatable :: theta_sfincs
nzpi = nz2pi/2
allocate (B_local(-nzpi:nzpi))
allocate (dBdz_local(-nzpi:nzpi))
allocate (gradpar_local(-nzpi:nzpi))
allocate (theta_sfincs(ntheta))
allocate (zed_stella(-nzpi:nzpi))
call init_zero_arrays
! first get some geometric quantities at this radius
! for theta from -pi to pi
q_local = geo_surf%qinp*(1.0+delrho*geo_surf%shat/geo_surf%rhoc)
B_local = bmag(1,-nzpi:nzpi) + delrho*dBdrho(-nzpi:nzpi)
dBdz_local = dbdzed(1,-nzpi:nzpi) + delrho*d2Bdrdth(-nzpi:nzpi)
gradpar_local = gradpar(-nzpi:nzpi) + delrho*dgradpardrho(-nzpi:nzpi)
zed_stella = zed(-nzpi:nzpi)+pi
theta_sfincs = export_f_theta(:ntheta)
iota = 1./q_local
! interpolate from stella zed-grid to sfincs theta grid
! point at -pi (stella) is same as point at 0 (sfincs)
BHat(1,1) = B_local(-nzpi)
call linear_interp_periodic (zed_stella, B_local, theta_sfincs(2:), BHat(2:,1))
! FLAG -- needs to be changed for stellarator runs
BHat = spread(BHat(:,1),2,nzeta)
dBHatdtheta(1,1) = dBdz_local(-nzpi)
call linear_interp_periodic (zed_stella, dBdz_local, theta_sfincs(2:), dBHatdtheta(2:,1))
dBHatdtheta = spread(dBHatdtheta(:,1),2,nzeta)
! this is bhat . grad theta
BHat_sup_theta(1,1) = B_local(-nzpi)*gradpar_local(-nzpi)
call linear_interp_periodic (zed_stella, B_local*gradpar_local, theta_sfincs(2:), BHat_sup_theta(2:,1))
BHat_sup_theta = spread(BHat_sup_theta(:,1),2,nzeta)
! this is I(psi) / (aref*Bref)
BHat_sub_zeta = geo_surf%rgeo + delrho*dIdrho
! this is grad psitor . (grad theta x grad zeta)
! note that + sign below relies on B = I grad zeta + grad zeta x grad psi
DHat = q_local*BHat_sup_theta
deallocate (B_local, dBdz_local, gradpar_local)
deallocate (theta_sfincs, zed_stella)
end subroutine pass_geometry_to_sfincs
subroutine init_zero_arrays
use globalVariables, only: dBHatdzeta
use globalVariables, only: dBHatdpsiHat
use globalVariables, only: BHat_sup_zeta
use globalVariables, only: BHat_sub_psi
use globalVariables, only: BHat_sub_theta
use globalVariables, only: dBHat_sub_psi_dtheta
use globalVariables, only: dBHat_sub_psi_dzeta
use globalVariables, only: dBHat_sub_theta_dpsiHat
use globalVariables, only: dBHat_sub_theta_dzeta
use globalVariables, only: dBHat_sub_zeta_dpsiHat
use globalVariables, only: dBHat_sub_zeta_dtheta
use globalVariables, only: dBHat_sup_theta_dpsiHat
use globalVariables, only: dBHat_sup_theta_dzeta
use globalVariables, only: dBHat_sup_zeta_dpsiHat
use globalVariables, only: dBHat_sup_zeta_dtheta
implicit none
dBHatdzeta = 0.
dBHatdpsiHat = 0.
BHat_sup_zeta = 0.
BHat_sub_psi = 0.
BHat_sub_theta = 0.
dBHat_sub_psi_dtheta = 0.
dBHat_sub_psi_dzeta = 0.
dBHat_sub_theta_dpsiHat = 0.
dBHat_sub_theta_dzeta = 0.
dBHat_sub_zeta_dpsiHat = 0.
dBHat_sub_zeta_dtheta = 0.
dBHat_sup_theta_dpsiHat = 0.
dBHat_sup_theta_dzeta = 0.
dBHat_sup_zeta_dpsiHat = 0.
dBHat_sup_zeta_dtheta = 0.
end subroutine init_zero_arrays
subroutine get_sfincs_output (f_neoclassical, phi_neoclassical, &
dfneo_dalpha, dphineo_dalpha)
use constants, only: pi
use sort, only: sort_array_ascending, unsort_array_ascending
use species, only: nspec
use zgrid, only: nzgrid, nz2pi
use export_f, only: h_sfincs => delta_f
use globalVariables, only: Phi1Hat
use kt_grids, only: nalpha
implicit none
real, dimension (:,-nzgrid:,:,:,:), intent (out) :: f_neoclassical
real, dimension (:,-nzgrid:), intent (out) :: phi_neoclassical
real, dimension (:,-nzgrid:,:,:,:), intent (out), optional :: dfneo_dalpha
real, dimension (:,-nzgrid:), intent (out), optional :: dphineo_dalpha
integer :: i, j
integer :: ialpha, iz, is
real, dimension (:), allocatable :: zed_stella
integer :: nfp, nfp_stella
integer, dimension (:), allocatable :: nzed_per_field_period
real, dimension (:,:), allocatable :: zed_stella_by_field_period
real, dimension (:,:), allocatable :: alpha_like_stella
integer, dimension (:,:), allocatable :: alpha_sort_map
integer :: nzed_sfincs, nalpha_sfincs
real, dimension (:), allocatable :: zed_sfincs, alpha_like_sfincs
real, dimension (:,:), allocatable :: phi_sfincs
real, dimension (:,:), allocatable :: phi_stella_zgrid, phi_stella
real, dimension (:,:), allocatable :: dphi_dalpha_stella_zgrid
real, dimension (:,:), allocatable :: tmp_sfincs, tmp_stella_zgrid, tmp_stella
real, dimension (:,:), allocatable :: dh_stella_zgrid
real, dimension (:,:,:,:), allocatable :: h_stella, dh_stella
! zed coordinate in stella is zeta when simulating stellarators (using vmec)
! and theta otherwise. this leads to some complications, treated below
allocate (zed_stella(nz2pi))
allocate (alpha_like_stella(nalpha,nz2pi))
allocate (alpha_sort_map(nalpha,nz2pi))
! obtain theta and zeta grids used in stella, and assign them
! to the zed and alpha-like coordinates.
! for stellarator, zed=zeta and alpha_like=theta.
! for tokamak, zed=theta and alpha_like = zeta.
call get_stella_theta_zeta_grids (alpha_like_stella, zed_stella)
! do iz = 1, nz2pi
! do ialpha = 1, size(alpha_like_stella,1)
! write (*,*) 'unsorted_alpha_like_stella', alpha_like_stella(ialpha,iz)
! end do
! write (*,*)
! end do
! rearrange alpha_like_stella to be in ascending order and store map
! so that sorting can be undone later
do iz = 1, nz2pi
call sort_array_ascending (alpha_like_stella(:,iz), alpha_sort_map(:,iz))
end do
! do iz = 1, nz2pi
! do ialpha = 1, size(alpha_like_stella,1)
! write (*,*) 'sorted_alpha_like_stella', alpha_like_stella(ialpha,iz)
! end do
! write (*,*)
! end do
! obtain the number of alpha-like and zed grid points to use in sfincs theta-zeta grid
! also obtain the number of field periods per 2*pi segment in zed
call get_sfincs_theta_zeta_grid_sizes (nalpha_sfincs, nzed_sfincs, nfp)
! obtain the alpha-like and zed coordinate grids
! note that additional points are added at periodic points
! that are not sampled in sfincs
allocate (zed_sfincs(nzed_sfincs))
allocate (alpha_like_sfincs(nalpha_sfincs))
call get_sfincs_theta_zeta_grids (alpha_like_sfincs, zed_sfincs)
allocate (phi_sfincs(nalpha_sfincs,nzed_sfincs))
call get_sfincs_field_theta_zeta (Phi1Hat, phi_sfincs)
! this is the number of field periods included in stella
! simulation domain
nfp_stella = int((zed_stella(nz2pi)*nfp-100.*epsilon(0.))/(2.*pi)) + 1
! obtain the number of zed grid points in each field period within the zed domain
allocate (nzed_per_field_period(nfp_stella))
call get_nzed_per_field_period (zed_stella, nfp, nzed_per_field_period)
! obtain the zed grid within each field period
allocate (zed_stella_by_field_period(maxval(nzed_per_field_period),nfp_stella))
call sort_zed_by_field_period (zed_stella, nzed_per_field_period, nfp, zed_stella_by_field_period)
! interpolate phi from sfincs zed grid to stella zed grid
allocate (phi_stella_zgrid(nalpha_sfincs,nz2pi))
call get_field_on_stella_zed_grid (phi_sfincs, nfp_stella, nfp, nzed_per_field_period, &
nalpha_sfincs, zed_sfincs, zed_stella_by_field_period, phi_stella_zgrid)
allocate (phi_stella(nalpha,nz2pi))
! interpolate onto stella (sorted) alpha grid
call get_field_stella (phi_stella_zgrid, alpha_like_sfincs, alpha_like_stella, phi_stella)
! need to remap from ascending (sorted) alpha grid to original ordering
do iz = 1, nz2pi
call unsort_array_ascending (phi_stella(:,iz), alpha_sort_map(:,iz))
end do
call get_field_on_extended_zed (phi_stella, phi_neoclassical)
if (present(dphineo_dalpha)) then
allocate (dphi_dalpha_stella_zgrid(nalpha_sfincs,nz2pi))
call get_dfield_dalpha (phi_stella_zgrid, alpha_like_sfincs, dphi_dalpha_stella_zgrid)
call get_field_stella (dphi_dalpha_stella_zgrid, alpha_like_sfincs, alpha_like_stella, phi_stella)
do iz = 1, nz2pi
call unsort_array_ascending (phi_stella(:,iz), alpha_sort_map(:,iz))
end do
call get_field_on_extended_zed (phi_stella, dphineo_dalpha)
deallocate (dphi_dalpha_stella_zgrid)
end if
deallocate (phi_stella, phi_stella_zgrid, phi_sfincs)
allocate (tmp_sfincs(nalpha_sfincs,nzed_sfincs))
allocate (tmp_stella_zgrid(nalpha_sfincs,nz2pi))
allocate (tmp_stella(nalpha,nz2pi))
allocate (h_stella(nalpha,nz2pi,size(h_sfincs,4),size(h_sfincs,5)))
if (present(dfneo_dalpha)) then
allocate (dh_stella_zgrid(nalpha_sfincs,nz2pi))
allocate (dh_stella(nalpha,nz2pi,size(h_sfincs,4),size(h_sfincs,5)))
end if
do is = 1, nspec
do i = 1, size(h_sfincs,5)
do j = 1, size(h_sfincs,4)
! re-order theta and zeta indices for sfincs h to ensure alpha-like coordinate
! appears before zed coordinate
call get_sfincs_field_theta_zeta (h_sfincs(is,:,:,j,i), tmp_sfincs)
! interpolate onto stella zed grid
call get_field_on_stella_zed_grid (tmp_sfincs, nfp_stella, nfp, nzed_per_field_period, &
nalpha_sfincs, zed_sfincs, zed_stella_by_field_period, tmp_stella_zgrid)
! interpolate onto (sorted) stella alpha-like coordinate
call get_field_stella (tmp_stella_zgrid, alpha_like_sfincs, alpha_like_stella, tmp_stella)
do iz = 1, nz2pi
call unsort_array_ascending (tmp_stella(:,iz), alpha_sort_map(:,iz))
end do
! use periodicity to copy onto extended zed grid if nperiod > 1
call get_field_on_extended_zed (tmp_stella, h_stella(:,:,j,i))
if (present(dfneo_dalpha)) then
call get_dfield_dalpha (tmp_stella_zgrid, alpha_like_sfincs, dh_stella_zgrid)
call get_field_stella (dh_stella_zgrid, alpha_like_sfincs, &
alpha_like_stella, tmp_stella)
do iz = 1, nz2pi
call unsort_array_ascending (tmp_stella(:,iz), alpha_sort_map(:,iz))
end do
call get_field_on_extended_zed (tmp_stella, dh_stella(:,:,j,i))
end if
end do
end do
do ialpha = 1, nalpha
do iz = -nzgrid, nzgrid
call sfincs_vspace_to_stella_vspace (ialpha, iz, is, h_stella(ialpha,iz+nzgrid+1,:,:), &
phi_neoclassical(ialpha,iz), f_neoclassical(ialpha,iz,:,:,is))
if (present(dfneo_dalpha)) call sfincs_vspace_to_stella_vspace (ialpha, iz, is, &
dh_stella(ialpha,iz+nzgrid+1,:,:), dphineo_dalpha(ialpha,iz), &
dfneo_dalpha(ialpha,iz,:,:,is))
end do
end do
end do
deallocate (tmp_sfincs, tmp_stella_zgrid, tmp_stella)
deallocate (zed_stella, alpha_like_stella)
deallocate (alpha_sort_map)
deallocate (zed_sfincs, alpha_like_sfincs)
deallocate (nzed_per_field_period, zed_stella_by_field_period)
deallocate (h_stella)
if (present(dfneo_dalpha)) deallocate (dh_stella, dh_stella_zgrid)
end subroutine get_sfincs_output
subroutine get_stella_theta_zeta_grids (alpha_like_stella, zed_stella)
use constants, only: pi
use zgrid, only: nz2pi, zed
use stella_geometry, only: zed_scalefac
use stella_geometry, only: alpha
use kt_grids, only: nalpha
use globalVariables, only: iota
implicit none
real, dimension (:,:), intent (out) :: alpha_like_stella
real, dimension (:), intent (out) :: zed_stella
integer :: nzpi
nzpi = nz2pi/2
! convert from scaled zed grid on [-pi,pi]