Skip to content

Commit 8579342

Browse files
authored
TokaMaker: Fix bugs with negative F0 and mode=1, fixes OpenFUSIONToolkit#321 (OpenFUSIONToolkit#322)
- Avoid new lint rules till handled - Add regression tests for negative F0 The error and fix effects the following subroutines and any downstream callers: - `grad_shaf_util.F90:save_ifile` - `grad_shaf_util.F90:save_eqdsk` - `grad_shaf_util.F90:sauter_fc` - `grad_shaf.F90:gs_get_chi` - `grad_shaf.F90:gs_save_fields` - `grad_shaf.F90:gs_get_qprof` - `grad_shaf.F90:gs_prof_interp_apply` - `grad_shaf.F90:gs_b_interp_apply` - `grad_shaf.F90:gs_save_prof` - `grad_shaf_fit.F90:fit_field_eval` - `tokamaker_f.F90:tokamaker_get_profs`
1 parent 3db6ce3 commit 8579342

6 files changed

Lines changed: 53 additions & 43 deletions

File tree

.github/workflows/lint.yaml

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -19,7 +19,7 @@ jobs:
1919
python3 -m venv ${{ github.workspace }}/oft_venv
2020
echo "source ${{ github.workspace }}/oft_venv/bin/activate" > ${{ github.workspace }}/setup_env.sh
2121
source ${{ github.workspace }}/setup_env.sh
22-
pip install ruff
22+
pip install "ruff<0.16.0"
2323
2424
- name: Check stack entries
2525
working-directory: src

src/physics/grad_shaf.F90

Lines changed: 6 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -3451,7 +3451,7 @@ subroutine gs_get_chi(self)
34513451
IF(self%mode==0)THEN
34523452
f=self%ffp_scale*self%I%f(psitmp(1))+self%I%f_offset
34533453
ELSE
3454-
f=SQRT(self%ffp_scale*self%I%F(psitmp(1)) + self%I%f_offset**2)
3454+
f=SIGN(1.d0,self%I%f_offset)*SQRT(self%ffp_scale*self%I%F(psitmp(1)) + self%I%f_offset**2)
34553455
END IF
34563456
do l=1,device%fe_rep%nce
34573457
call oft_blag_eval(device%fe_rep,j,l,device%fe_rep%quad%pts(:,m),rop)
@@ -4084,7 +4084,7 @@ subroutine gs_save_fields(self,pts,npts,filename)
40844084
IF(self%mode==0)THEN
40854085
B(2)=(self%ffp_scale*self%I%f(psitmp(1))+self%I%f_offset)/pts(1,i)
40864086
ELSE
4087-
B(2)=SQRT(self%ffp_scale*self%I%f(psitmp(1)) + self%I%f_offset**2)/pts(1,i)
4087+
B(2)=SIGN(1.d0,self%I%f_offset)*SQRT(self%ffp_scale*self%I%f(psitmp(1)) + self%I%f_offset**2)/pts(1,i)
40884088
END IF
40894089
B=B*self%psiscale
40904090
! Handle anisotropic pressure
@@ -4477,8 +4477,7 @@ subroutine gs_get_qprof(gseq,nr,psi_q,prof,dl,rbounds,zbounds,ravgs)
44774477
IF(gseq%mode==0)THEN
44784478
fpol=gseq%ffp_scale*gseq%I%f(psi_surf)+gseq%I%f_offset
44794479
ELSE
4480-
fpol=SQRT(gseq%ffp_scale*gseq%I%f(psi_surf) + gseq%I%f_offset**2) &
4481-
+ gseq%I%f_offset*(1.d0-SIGN(1.d0,gseq%I%f_offset))
4480+
fpol=SIGN(1.d0,gseq%I%f_offset)*SQRT(gseq%ffp_scale*gseq%I%f(psi_surf) + gseq%I%f_offset**2)
44824481
END IF
44834482
!---Safety Factor (q)
44844483
qpsi=fpol*active_tracer%v(3)/(2*pi)
@@ -4640,7 +4639,7 @@ subroutine gs_prof_interp_apply(self,cell,f,gop,val)
46404639
IF(self%equil%mode==0)THEN
46414640
val(1)=self%equil%psiscale*(self%equil%ffp_scale*self%equil%I%f(psitmp(1))+self%equil%I%f_offset)
46424641
ELSE
4643-
val(1)=self%equil%psiscale*SQRT(self%equil%ffp_scale*self%equil%I%f(psitmp(1)) + self%equil%I%f_offset**2)
4642+
val(1)=self%equil%psiscale*SIGN(1.d0,self%equil%I%f_offset)*SQRT(self%equil%ffp_scale*self%equil%I%f(psitmp(1)) + self%equil%I%f_offset**2)
46444643
END IF
46454644
ELSE
46464645
val(1)=self%equil%psiscale*self%equil%I%f_offset
@@ -4698,7 +4697,7 @@ subroutine gs_b_interp_apply(self,cell,f,gop,val)
46984697
IF(self%equil%mode==0)THEN
46994698
val(2)=self%equil%psiscale*(self%equil%ffp_scale*self%equil%I%f(psitmp(1))+self%equil%I%f_offset)/(pt(1)+gs_epsilon)
47004699
ELSE
4701-
val(2)=self%equil%psiscale*SQRT(self%equil%ffp_scale*self%equil%I%f(psitmp(1)) + self%equil%I%f_offset**2)/(pt(1)+gs_epsilon)
4700+
val(2)=self%equil%psiscale*SIGN(1.d0,self%equil%I%f_offset)*SQRT(self%equil%ffp_scale*self%equil%I%f(psitmp(1)) + self%equil%I%f_offset**2)/(pt(1)+gs_epsilon)
47024701
END IF
47034702
ELSE
47044703
val(2)=self%equil%psiscale*self%equil%I%f_offset/(pt(1)+gs_epsilon)
@@ -5156,7 +5155,7 @@ subroutine gs_save_prof(self,filename,mpsi_sample)
51565155
outtmp(2)=self%ffp_scale*self%I%fp(r)
51575156
outtmp(3)=self%psiscale*self%ffp_scale*self%I%f(r) + self%I%f_offset
51585157
ELSE
5159-
outtmp(3)=SQRT(self%psiscale*self%ffp_scale*self%I%f(r) + self%I%f_offset**2)
5158+
outtmp(3)=SIGN(1.d0,self%I%f_offset)*SQRT(self%psiscale*self%ffp_scale*self%I%f(r) + self%I%f_offset**2)
51605159
outtmp(2)=self%ffp_scale*self%I%fp(r)/(2.d0*outtmp(3))
51615160
END IF
51625161
outtmp(4)=self%psiscale*self%p_scale*self%P%fp(r)

src/physics/grad_shaf_fit.F90

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -537,20 +537,20 @@ function minpack_exit_reason(info) result(exit_reason)
537537
! is at most xtol.
538538
!
539539
! info = 3 conditions for info = 1 and info = 2 both hold.
540-
!
540+
!
541541
! info = 4 the cosine of the angle between fvec and any
542542
! column of the jacobian is at most gtol in
543543
! absolute value.
544-
!
544+
!
545545
! info = 5 number of calls to fcn with iflag = 1 has
546546
! reached maxfev.
547-
!
547+
!
548548
! info = 6 ftol is too small. no further reduction in
549549
! the sum of squares is possible.
550-
!
550+
!
551551
! info = 7 xtol is too small. no further improvement in
552552
! the approximate solution x is possible.
553-
!
553+
!
554554
! info = 8 gtol is too small. fvec is orthogonal to the
555555
! columns of the jacobian to machine precision.
556556
SELECT CASE(info)
@@ -1444,7 +1444,7 @@ FUNCTION fit_field_eval(self,gs) RESULT(val)
14441444
IF(gs%mode==0)THEN
14451445
btmp(2)= gs%psiscale*(gs%ffp_scale*gs%I%f(psi(1))+gs%I%f_offset)/self%r(1)
14461446
ELSE
1447-
btmp(2)=gs%psiscale*SQRT(gs%ffp_scale*gs%I%f(psi(1)) + gs%I%f_offset**2)/self%r(1)
1447+
btmp(2)=gs%psiscale*SIGN(1.d0,gs%I%f_offset)*SQRT(gs%ffp_scale*gs%I%f(psi(1)) + gs%I%f_offset**2)/self%r(1)
14481448
END IF
14491449
btmp(3)= gs%psiscale*gpsi(1)/self%r(1)
14501450
!---

src/physics/grad_shaf_util.F90

Lines changed: 5 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -867,8 +867,7 @@ subroutine gs_save_ifile(gseq,filename,npsi,ntheta,psi_pad,lcfs_press,pack_lcfs,
867867
IF(gseq%mode==0)THEN
868868
cout(j,2)=gseq%ffp_scale*gseq%I%f(psi_surf(1))+gseq%I%f_offset
869869
ELSE
870-
cout(j,2)=SQRT(gseq%ffp_scale*gseq%I%f(psi_surf(1)) + gseq%I%f_offset**2) &
871-
+ gseq%I%f_offset*(1.d0-SIGN(1.d0,gseq%I%f_offset))
870+
cout(j,2)=SIGN(1.d0,gseq%I%f_offset)*SQRT(gseq%ffp_scale*gseq%I%f(psi_surf(1)) + gseq%I%f_offset**2)
872871
END IF
873872
cout(j,3)=gseq%p_scale*gseq%P%f(psi_surf(1))/mu0 ! Plasma pressure
874873
cout(j,4)=cout(j,2)*active_tracer%v(3)/(2*pi) ! Safety Factor (q)
@@ -891,8 +890,7 @@ subroutine gs_save_ifile(gseq,filename,npsi,ntheta,psi_pad,lcfs_press,pack_lcfs,
891890
IF(gseq%mode==0)THEN
892891
cout(1,2)=(gseq%ffp_scale*gseq%I%f(x2)+gseq%I%f_offset)
893892
ELSE
894-
cout(1,2)=SQRT(gseq%ffp_scale*gseq%I%f(x2) + gseq%I%f_offset**2) &
895-
+ gseq%I%f_offset*(1.d0-SIGN(1.d0,gseq%I%f_offset))
893+
cout(1,2)=SIGN(1.d0,gseq%I%f_offset)*SQRT(gseq%ffp_scale*gseq%I%f(x2) + gseq%I%f_offset**2)
896894
END IF
897895
cout(1,3)=gseq%p_scale*gseq%P%f(x2)/mu0
898896
cout(1,4)=(cout(3,4)-cout(2,4))*(x2-cout(2,1))/(cout(3,1)-cout(2,1)) + cout(2,4)
@@ -1132,10 +1130,8 @@ subroutine gs_save_eqdsk(gseq,filename,nr,nz,rbounds,zbounds,run_info,limiter_fi
11321130
fpol(j)=gseq%ffp_scale*gseq%I%f(psi_surf)+gseq%I%f_offset
11331131
ffprim(j)=gseq%I%fp(psi_surf)*((gseq%ffp_scale**2)*gseq%I%f(psi_surf)+gseq%ffp_scale*gseq%I%f_offset)
11341132
ELSE
1135-
fptmp=SQRT(gseq%ffp_scale*gseq%I%f(psi_trace) + gseq%I%f_offset**2) &
1136-
+ gseq%I%f_offset*(1.d0-SIGN(1.d0,gseq%I%f_offset))
1137-
fpol(j)=SQRT(gseq%ffp_scale*gseq%I%f(psi_surf) + gseq%I%f_offset**2) &
1138-
+ gseq%I%f_offset*(1.d0-SIGN(1.d0,gseq%I%f_offset))
1133+
fptmp=SIGN(1.d0,gseq%I%f_offset)*SQRT(gseq%ffp_scale*gseq%I%f(psi_trace) + gseq%I%f_offset**2)
1134+
fpol(j)=SIGN(1.d0,gseq%I%f_offset)*SQRT(gseq%ffp_scale*gseq%I%f(psi_surf) + gseq%I%f_offset**2)
11391135
ffprim(j)=0.5d0*gseq%ffp_scale*gseq%I%fp(psi_surf)
11401136
END IF
11411137
pres(j)=gseq%p_scale*gseq%P%f(psi_surf)/mu0
@@ -1414,8 +1410,7 @@ subroutine sauter_fc(gseq,nr,psi_q,fc,r_avgs,modb_avgs)
14141410
IF(gseq%mode==0)THEN
14151411
field%f_surf=gseq%ffp_scale*gseq%I%f(psi_surf)+gseq%I%f_offset
14161412
ELSE
1417-
field%f_surf=SQRT(gseq%ffp_scale*gseq%I%f(psi_surf) + gseq%I%f_offset**2) &
1418-
+ gseq%I%f_offset*(1.d0-SIGN(1.d0,gseq%I%f_offset))
1413+
field%f_surf=SIGN(1.d0,gseq%I%f_offset)*SQRT(gseq%ffp_scale*gseq%I%f(psi_surf) + gseq%I%f_offset**2)
14191414
END IF
14201415
field%bmax=0.d0
14211416
field%stage_1=.TRUE.

src/python/wrappers/tokamaker_f.F90

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1302,7 +1302,7 @@ SUBROUTINE tokamaker_get_profs(tMaker_equil_ptr,npsi,psi_in,f,fp,p,pp,error_str)
13021302
fp(i)=tMaker_equil_obj%ffp_scale*tMaker_equil_obj%I%fp(r)
13031303
f(i)=tMaker_equil_obj%psiscale*tMaker_equil_obj%ffp_scale*tMaker_equil_obj%I%f(r) + tMaker_equil_obj%I%f_offset
13041304
ELSE
1305-
f(i)=SQRT(tMaker_equil_obj%psiscale*tMaker_equil_obj%ffp_scale*tMaker_equil_obj%I%f(r) + tMaker_equil_obj%I%f_offset**2)
1305+
f(i)=SIGN(1.d0,tMaker_equil_obj%I%f_offset)*SQRT(tMaker_equil_obj%psiscale*tMaker_equil_obj%ffp_scale*tMaker_equil_obj%I%f(r) + tMaker_equil_obj%I%f_offset**2)
13061306
fp(i)=tMaker_equil_obj%ffp_scale*tMaker_equil_obj%I%fp(r)/(2.d0*f(i))
13071307
END IF
13081308
pp(i)=tMaker_equil_obj%psiscale*tMaker_equil_obj%p_scale*tMaker_equil_obj%P%fp(r)

src/tests/physics/test_TokaMaker.py

Lines changed: 34 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -86,7 +86,7 @@ def validate_dict(results,dict_exp,tol_dict=None):
8686
return test_result
8787

8888

89-
def validate_eqdsk(file_test,file_ref):
89+
def validate_eqdsk(file_test,file_ref,helicity=1.0):
9090
from OpenFUSIONToolkit.TokaMaker.util import read_eqdsk
9191
try:
9292
test_data = read_eqdsk(file_test)
@@ -98,6 +98,9 @@ def validate_eqdsk(file_test,file_ref):
9898
except:
9999
print("FAILED: Could not read reference EQDSK")
100100
return False
101+
ref_data['bcentr'] *= helicity
102+
ref_data['fpol'] *= helicity
103+
ref_data['qpsi'] *= helicity
101104
test_result = True
102105
for key, exp_val in ref_data.items():
103106
result_val = test_data.get(key,None)
@@ -121,7 +124,7 @@ def validate_eqdsk(file_test,file_ref):
121124
return test_result
122125

123126

124-
def validate_ifile(ifile_test,ifile_ref):
127+
def validate_ifile(ifile_test,ifile_ref,helicity=1.0):
125128
from OpenFUSIONToolkit.TokaMaker.util import read_ifile
126129
try:
127130
test_data = read_ifile(ifile_test)
@@ -133,6 +136,8 @@ def validate_ifile(ifile_test,ifile_ref):
133136
except:
134137
print("FAILED: Could not read reference i-file")
135138
return False
139+
ref_data['f'] *= helicity
140+
ref_data['q'] *= helicity
136141
test_result = True
137142
for key, exp_val in ref_data.items():
138143
result_val = test_data.get(key,None)
@@ -467,7 +472,7 @@ def test_coil_h3(order,dist_coil):
467472

468473

469474
#============================================================================
470-
def run_ITER_case(mesh_resolution,fe_orders,test_type,mp_q):
475+
def run_ITER_case(mesh_resolution,fe_orders,test_type,helicity,mp_q):
471476
def create_mesh():
472477
with open('ITER_geom.json','r') as fid:
473478
ITER_geom = json.load(fid)
@@ -515,7 +520,7 @@ def create_mesh():
515520
mesh_pts,mesh_lc,mesh_reg,coil_dict,cond_dict = load_gs_mesh('ITER_mesh.h5')
516521
mygs.setup_mesh(mesh_pts,mesh_lc,mesh_reg)
517522
mygs.setup_regions(cond_dict=cond_dict,coil_dict=coil_dict)
518-
mygs.setup(order=fe_order,F0=5.3*6.2)
523+
mygs.setup(order=fe_order,F0=helicity*5.3*6.2)
519524
#
520525
if test_type.startswith('eig'):
521526
if test_type == 'eig_dep':
@@ -740,10 +745,10 @@ def test_ITER_eig(order):
740745
exp_dict = {
741746
'Tau_w': [6.619977E-01, 3.479492E-01, 2.554444E-01, 1.910381E-01, 1.782464E-01]
742747
}
743-
results = mp_run(run_ITER_case,(1.0,(order,),'eig'))
748+
results = mp_run(run_ITER_case,(1.0,(order,),'eig',1.0))
744749
assert validate_dict(results,exp_dict)
745750
# Test deprecated interface
746-
results = mp_run(run_ITER_case,(1.0,(order,),'eig_dep'))
751+
results = mp_run(run_ITER_case,(1.0,(order,),'eig_dep',1.0))
747752
assert validate_dict(results,exp_dict)
748753

749754
@pytest.mark.coverage
@@ -753,10 +758,10 @@ def test_ITER_stability(order):
753758
'gamma': [12.3620, -1.83981, -3.41613, -5.12470, -6.53393],
754759
'nl_change': [225.4421413167051, 338.0113029638385][order-2]
755760
}
756-
results = mp_run(run_ITER_case,(1.0,(order,),'stab'))
761+
results = mp_run(run_ITER_case,(1.0,(order,),'stab',1.0))
757762
assert validate_dict(results,exp_dict)
758763
# Test deprecated interface
759-
results = mp_run(run_ITER_case,(1.0,(order,),'stab_dep'))
764+
results = mp_run(run_ITER_case,(1.0,(order,),'stab_dep',1.0))
760765
assert validate_dict(results,exp_dict)
761766

762767
ITER_eq_dict = {
@@ -788,18 +793,29 @@ def test_ITER_stability(order):
788793

789794
@pytest.mark.coverage
790795
@pytest.mark.parametrize("order", (2,3))#,4))
791-
def test_ITER_eq(order):
792-
results = mp_run(run_ITER_case,(1.0,(order,),''))
793-
assert validate_dict(results,ITER_eq_dict)
794-
assert validate_eqdsk('tokamaker.eqdsk','ITER_test.eqdsk')
795-
assert validate_ifile('tokamaker.ifile','ITER_test.ifile')
796+
@pytest.mark.parametrize("helicity", (1.0,-1.0))
797+
def test_ITER_eq(order,helicity):
798+
eq_dict = ITER_eq_dict.copy()
799+
eq_dict['tflux'] *= helicity
800+
eq_dict['dflux'] *= helicity
801+
eq_dict['q_0'] *= helicity
802+
eq_dict['q_95'] *= helicity
803+
results = mp_run(run_ITER_case,(1.0,(order,),'',helicity))
804+
assert validate_dict(results,eq_dict)
805+
assert validate_eqdsk('tokamaker.eqdsk','ITER_test.eqdsk',helicity)
806+
assert validate_ifile('tokamaker.ifile','ITER_test.ifile',helicity)
796807

797808
@pytest.mark.coverage
798809
@pytest.mark.parametrize("order", (2,3))#,4))
799-
def test_ITER_eq_io(order):
810+
@pytest.mark.parametrize("helicity", (1.0,-1.0))
811+
def test_ITER_eq_io(order,helicity):
800812
eq_dict = ITER_eq_dict.copy()
801813
eq_dict['nl_its'] = 1
802-
results = mp_run(run_ITER_case,(1.0,(order,),'io'))
814+
eq_dict['tflux'] *= helicity
815+
eq_dict['dflux'] *= helicity
816+
eq_dict['q_0'] *= helicity
817+
eq_dict['q_95'] *= helicity
818+
results = mp_run(run_ITER_case,(1.0,(order,),'io',helicity))
803819
assert validate_dict(results,eq_dict)
804820

805821
@pytest.mark.coverage
@@ -813,7 +829,7 @@ def test_ITER_recon():
813829
ITER_recon_dict['l_i'] = 0.8872565
814830
ITER_recon_dict['beta_tor'] = 1.906180
815831
ITER_recon_dict['beta_n'] = 1.284312
816-
results = mp_run(run_ITER_case,(1.0,(2,),'recon'))
832+
results = mp_run(run_ITER_case,(1.0,(2,),'recon',1.0))
817833
assert validate_dict(results,ITER_recon_dict)
818834

819835
@pytest.mark.coverage
@@ -829,11 +845,11 @@ def test_ITER_recon_legacy():
829845
ITER_recon_dict['l_i'] = 0.8906732
830846
ITER_recon_dict['beta_tor'] = 1.934024
831847
ITER_recon_dict['beta_n'] = 1.294385
832-
results = mp_run(run_ITER_case,(1.0,(2,),'recon_legacy'))
848+
results = mp_run(run_ITER_case,(1.0,(2,),'recon_legacy',1.0))
833849
assert validate_dict(results,ITER_recon_dict)
834850

835851
def test_ITER_concurrent():
836-
results = mp_run(run_ITER_case,(1.0,(2,3),''))
852+
results = mp_run(run_ITER_case,(1.0,(2,3),'',1.0))
837853
assert validate_dict(results,ITER_eq_dict)
838854

839855
#============================================================================

0 commit comments

Comments
 (0)