Skip to content

Commit a8836f2

Browse files
author
Caoxiang Zhu
committed
some minor changes
1 parent 945f95e commit a8836f2

3 files changed

Lines changed: 213 additions & 13 deletions

File tree

‎New/macros‎

Lines changed: 38 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -292,6 +292,44 @@ m4_define(HWRITERC,{
292292

293293
})m4_dnl
294294

295+
m4_define(HWRITECH,{
296+
297+
!macro expansion of hwritech; write character in hdf5 format;
298+
299+
rank = 1 ; onedims(1) = $1
300+
301+
if( $1.le.0 ) then
302+
303+
write(0,'("macros : fatal : error calling hwriterv ; $2 : $1.le.0 ;")')
304+
305+
else
306+
307+
call h5screate_simple_f( rank , onedims , space_id , hdfier )
308+
if( hdfier.ne.0 ) then
309+
write(0,'("macros : fatal : error calling h5screate_simple_f ;")')
310+
call MPI_ABORT( MPI_COMM_WORLD, 1, ierr )
311+
stop "macros : error calling h5screate_simple_f;"
312+
endif
313+
314+
call h5dcreate_f( file_id , "$2" , H5T_NATIVE_CHARACTER , space_id , dset_id , hdfier )
315+
if( hdfier.ne.0 ) then
316+
write(0,'("macros : fatal : error calling h5dcreate_f ;")')
317+
call MPI_ABORT( MPI_COMM_WORLD, 1, ierr )
318+
stop "macros : error calling h5dcreate_f;"
319+
endif
320+
321+
call h5dwrite_f( dset_id , H5T_NATIVE_CHARACTER , $3 , onedims , hdfier )
322+
if( hdfier.ne.0 ) then
323+
write(0,'("macros : fatal : error calling h5dwrite_f ;")')
324+
call MPI_ABORT( MPI_COMM_WORLD, 1, ierr )
325+
stop "macros : error calling h5dwrite_f;"
326+
endif
327+
328+
CALL h5dclose_f(dset_id, hdfier) ! terminate dataset;
329+
330+
endif
331+
332+
})m4_dnl
295333

296334
m4_define(TMPOUT,{ !Temperarily output a message to help debugging
297335
if( myid .eq. 0) write(ounit,*) "### DEBUG : $1 = ", $1

‎New/saving.h‎

Lines changed: 5 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -29,7 +29,8 @@ subroutine saving
2929

3030

3131
INTEGER :: ii, jj, icoil, NF
32-
! REAL, allocatable :: perA(:)
32+
CHARACTER(LEN=10) :: version='v0.1.0'
33+
3334

3435
! the following are used by the macros HWRITEXX below; do not alter/remove;
3536
INTEGER :: hdfier, rank
@@ -65,6 +66,8 @@ subroutine saving
6566
call h5fcreate_f( trim(hdf5file), H5F_ACC_TRUNC_F, file_id, hdfier ) ! create new file;
6667
FATAL( restart, hdfier.ne.0, error calling h5fcreate_f )
6768

69+
HWRITECH( 10 , version , version )
70+
6871
!INPUT namelist;
6972
HWRITEIV( 1 , IsQuiet , IsQuiet )
7073
HWRITEIV( 1 , IsSymmetric , IsSymmetric )
@@ -82,6 +85,7 @@ subroutine saving
8285
HWRITEIV( 1 , NFcoil , NFcoil )
8386
HWRITEIV( 1 , Nseg , Nseg )
8487
HWRITEIV( 1 , case_optimize , case_optimize )
88+
HWRITERV( 1 , exit_xtol , exit_tol )
8589
HWRITEIV( 1 , IsNormalize , IsNormalize )
8690
HWRITEIV( 1 , ISNormWeight , IsNormWeight )
8791
HWRITEIV( 1 , case_bnormal , case_bnormal )

‎pyfocus/coil.py‎

Lines changed: 170 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -19,6 +19,8 @@
1919
12. plot coils from FOCUS hdf5 output data;
2020
13. animating the coils evolution from FOCUS hdf5 output data;
2121
14. plot cost function converging curves from FOCUS hdf5 output data;
22+
'''
23+
2224
'''
2325
#QT_API imcompactible
2426
import sip
@@ -31,6 +33,7 @@
3133
from PyQt4.QtSvg import *
3234
from PyQt4.QtCore import pyqtSignal as Signal
3335
from PyQt4.QtCore import pyqtSlot as Slot
36+
'''
3437

3538
import numpy as np
3639
import matplotlib.pyplot as plt
@@ -66,6 +69,8 @@ def __init__(self, hdf5file):
6669
self.names = hdf5.keys()
6770
for name in self.names:
6871
setattr(self, name, hdf5[name].value)
72+
if hdf5[name].value.dtype == 'S1': #only for hdf5 string;
73+
setattr(self, name, ''.join(hdf5[name].value))
6974
except IOError:
7075
print "Error: this is not a valid hdf5 file."
7176
else:
@@ -656,7 +661,7 @@ def plot_poincare( path, pp, code='Glass', prange='full', dotsize=0.5):
656661

657662
#------------------------------------------ 12 ----------------------------------------------------
658663

659-
def hdf5coil(h5data, nn = -1):
664+
def hdf5coil(h5data, nn = -1, old=False):
660665
'''
661666
read FOCUS hdf5 output file and return the nn-th coil data
662667
'''
@@ -667,24 +672,42 @@ def hdf5coil(h5data, nn = -1):
667672
assert ncoil>0
668673
nf = h5data.NFcoil[0]
669674
assert nf>0
670-
nseg = h5data.Nseg[0]
675+
if old:
676+
nseg = h5data.NDcoil[0]
677+
else:
678+
nseg = h5data.Nseg[0]
671679
assert nseg>0
672680
coilfou = np.reshape(h5data.coilspace[:,nn], (ncoil,-1)) #all coil data at nn time; and reshape
673681
coildata = np.ndarray((ncoil,),dtype=np.object)
674682

675683
for icoil in range(ncoil):
676-
xyzfou = np.reshape(coilfou[icoil,1:], (3,-1)) #pop current and reshape in 3 rows;
684+
685+
if old: #old format hdf5
686+
tmpfou = np.reshape(coilfou[icoil, 4:], (nf, 6))
687+
xyzfou = np.zeros([3, 2*nf+1], dtype=np.float32)
688+
xyzfou[0, 0] = coilfou[icoil, 1] #xc0
689+
xyzfou[1, 0] = coilfou[icoil, 2] #yc0
690+
xyzfou[2, 0] = coilfou[icoil, 3] #yc0
691+
xyzfou[0, 1:nf+1] = tmpfou[:, 0] #xc(1:n)
692+
xyzfou[1, 1:nf+1] = tmpfou[:, 1] #yc(1:n)
693+
xyzfou[2, 1:nf+1] = tmpfou[:, 2] #zc(1:n)
694+
xyzfou[0, nf+1:2*nf+1] = tmpfou[:, 3] #xs(1:n)
695+
xyzfou[1, nf+1:2*nf+1] = tmpfou[:, 4] #ys(1:n)
696+
xyzfou[2, nf+1:2*nf+1] = tmpfou[:, 5] #zs(1:n)
697+
else:
698+
xyzfou = np.reshape(coilfou[icoil,1:], (3,-1)) #pop current and reshape in 3 rows;
699+
677700
coildata[icoil] = coil()
678701
#coildata[icoil].I = coilfou[icoil,0] #coil current
679702
coildata[icoil].x.extend(fourier(xyzfou[0,:]).disc(nseg))
680703
coildata[icoil].y.extend(fourier(xyzfou[1,:]).disc(nseg))
681704
coildata[icoil].z.extend(fourier(xyzfou[2,:]).disc(nseg))
682-
print "Read {} coils in {}.".format(ncoil, filename)
705+
print "Read {} coils.".format(ncoil)
683706
return coildata
684707

685708
#------------------------------------------ 13 ----------------------------------------------------
686709

687-
def coilevolve(h5data, delay = 250):
710+
def coilevolve(h5data, delay = 250, nout=-1, old = False, savepng=False):
688711
'''
689712
plot coil evolution movie
690713
'''
@@ -696,11 +719,21 @@ def coilevolve(h5data, delay = 250):
696719
assert ncoil>0
697720
nf = h5data.NFcoil[0]
698721
assert nf>0
699-
nseg = h5data.Nseg[0]
722+
if old:
723+
nseg = h5data.NDcoil[0]
724+
else:
725+
nseg = h5data.Nseg[0]
700726
assert nseg>0
701-
nout = h5data.iout[0]
727+
if nout < 0:
728+
if old:
729+
nout = h5data.itau[0]
730+
else:
731+
nout = h5data.iout[0]
702732
assert nout>1
703-
nfp = h5data.Nfp[0]
733+
if old:
734+
nfp = 1
735+
else:
736+
nfp = h5data.Nfp[0]
704737
unicoil = ncoil/nfp
705738

706739
xx = np.zeros([ncoil,nseg], dtype=np.float32)
@@ -712,34 +745,82 @@ def coilevolve(h5data, delay = 250):
712745
#initial coils (fixed for comparison)
713746
coilinit = np.reshape(h5data.coilspace[:,0], (ncoil,-1))
714747
for icoil in range(ncoil):
715-
xyzfou = np.reshape(coilinit[icoil,1:], (3,-1))
748+
if old: #old format hdf5
749+
tmpfou = np.reshape(coilinit[icoil, 4:], (nf, 6))
750+
xyzfou = np.zeros([3, 2*nf+1], dtype=np.float32)
751+
xyzfou[0, 0] = coilinit[icoil, 1] #xc0
752+
xyzfou[1, 0] = coilinit[icoil, 2] #yc0
753+
xyzfou[2, 0] = coilinit[icoil, 3] #yc0
754+
xyzfou[0, 1:nf+1] = tmpfou[:, 0] #xc(1:n)
755+
xyzfou[1, 1:nf+1] = tmpfou[:, 1] #yc(1:n)
756+
xyzfou[2, 1:nf+1] = tmpfou[:, 2] #zc(1:n)
757+
xyzfou[0, nf+1:2*nf+1] = tmpfou[:, 3] #xs(1:n)
758+
xyzfou[1, nf+1:2*nf+1] = tmpfou[:, 4] #ys(1:n)
759+
xyzfou[2, nf+1:2*nf+1] = tmpfou[:, 5] #zs(1:n)
760+
else:
761+
xyzfou = np.reshape(coilinit[icoil,1:], (3,-1)) #pop current and reshape in 3 rows;
762+
716763
mlab.plot3d(fourier(xyzfou[0,:]).disc(nseg), fourier(xyzfou[1,:]).disc(nseg),
717764
fourier(xyzfou[2,:]).disc(nseg), color=(0.5, 0.5, 0.5))
718765

719766
# first coils
720767
coilfou = np.reshape(h5data.coilspace[:,0], (ncoil,-1))
721768
for icoil in range(ncoil):
722-
xyzfou = np.reshape(coilfou[icoil,1:], (3,-1))
769+
if old: #old format hdf5
770+
tmpfou = np.reshape(coilfou[icoil, 4:], (nf, 6))
771+
xyzfou = np.zeros([3, 2*nf+1], dtype=np.float32)
772+
xyzfou[0, 0] = coilfou[icoil, 1] #xc0
773+
xyzfou[1, 0] = coilfou[icoil, 2] #yc0
774+
xyzfou[2, 0] = coilfou[icoil, 3] #yc0
775+
xyzfou[0, 1:nf+1] = tmpfou[:, 0] #xc(1:n)
776+
xyzfou[1, 1:nf+1] = tmpfou[:, 1] #yc(1:n)
777+
xyzfou[2, 1:nf+1] = tmpfou[:, 2] #zc(1:n)
778+
xyzfou[0, nf+1:2*nf+1] = tmpfou[:, 3] #xs(1:n)
779+
xyzfou[1, nf+1:2*nf+1] = tmpfou[:, 4] #ys(1:n)
780+
xyzfou[2, nf+1:2*nf+1] = tmpfou[:, 5] #zs(1:n)
781+
else:
782+
xyzfou = np.reshape(coilfou[icoil,1:], (3,-1)) #pop current and reshape in 3 rows;
783+
723784
xx[icoil,:] = fourier(xyzfou[0,:]).disc(nseg)
724785
yy[icoil,:] = fourier(xyzfou[1,:]).disc(nseg)
725786
zz[icoil,:] = fourier(xyzfou[2,:]).disc(nseg)
726787
l = mlab.plot3d(xx[icoil,:], yy[icoil,:], zz[icoil,:], color=c[icoil%unicoil])
727788
coils.append(l)
728789

729-
@mlab.show
790+
print ' in coilevolve : '
730791
@mlab.animate(delay=delay)
731792
def anim():
793+
print ' in anim '
732794
while 1:
733795
for iout in range(nout):
734796
coilfou = np.reshape(h5data.coilspace[:,iout], (ncoil,-1)) #all coil data at iout time;
735797
for icoil in range(ncoil):
736-
xyzfou = np.reshape(coilfou[icoil,1:], (3,-1))
798+
if old: #old format hdf5
799+
tmpfou = np.reshape(coilfou[icoil, 4:], (nf, 6))
800+
xyzfou = np.zeros([3, 2*nf+1], dtype=np.float32)
801+
xyzfou[0, 0] = coilfou[icoil, 1] #xc0
802+
xyzfou[1, 0] = coilfou[icoil, 2] #yc0
803+
xyzfou[2, 0] = coilfou[icoil, 3] #yc0
804+
xyzfou[0, 1:nf+1] = tmpfou[:, 0] #xc(1:n)
805+
xyzfou[1, 1:nf+1] = tmpfou[:, 1] #yc(1:n)
806+
xyzfou[2, 1:nf+1] = tmpfou[:, 2] #zc(1:n)
807+
xyzfou[0, nf+1:2*nf+1] = tmpfou[:, 3] #xs(1:n)
808+
xyzfou[1, nf+1:2*nf+1] = tmpfou[:, 4] #ys(1:n)
809+
xyzfou[2, nf+1:2*nf+1] = tmpfou[:, 5] #zs(1:n)
810+
else:
811+
xyzfou = np.reshape(coilfou[icoil,1:], (3,-1)) #pop current and reshape in 3 rows;
812+
737813
xx[icoil,:] = fourier(xyzfou[0,:]).disc(nseg)
738814
yy[icoil,:] = fourier(xyzfou[1,:]).disc(nseg)
739815
zz[icoil,:] = fourier(xyzfou[2,:]).disc(nseg)
740816
coils[icoil].mlab_source.set(x=xx[icoil,:], y=yy[icoil,:], z=zz[icoil,:])
817+
if savepng :
818+
filename = 'coil_evolution_'+str(iout).zfill(6)+'.png'
819+
print(filename)
820+
mlab.savefig(filename=filename)
741821
yield
742822
anim()
823+
mlab.show()
743824
return
744825

745826
#------------------------------------------ 14 ----------------------------------------------------
@@ -793,3 +874,80 @@ def chievolve(h5data, term='chi', width = 2.5, linestyle='-', color='b'):
793874
plt.legend(loc='upper right', frameon=False, prop={'size':24, 'weight':'bold'})
794875

795876
return
877+
878+
#------------------------------------------ 15 ----------------------------------------------------
879+
880+
def plot_hdf5_surface(h5data, plottype='surf', color=(1,0,0)):
881+
'''
882+
plot plasma surface from the hdf5 data
883+
884+
plottype:
885+
'surf' : just plot the surface with pure color;
886+
'targetBn' : plot the surface with target Bn distribution;
887+
'realizedBn': plot the surface with realized Bn distribution;
888+
'''
889+
890+
if not isinstance(h5data, hdf5):
891+
print h5data, " is not a valid FOCUS hdf5 file."
892+
return
893+
894+
# data manipulation
895+
nt = h5data.Nteta[0] #opposite notation for nt & nz!
896+
nz = h5data.Nzeta[0]
897+
xx = np.zeros([nt+1,nz+1], dtype=np.float32)
898+
yy = np.zeros([nt+1,nz+1], dtype=np.float32)
899+
zz = np.zeros([nt+1,nz+1], dtype=np.float32)
900+
Bn = np.zeros([nt+1,nz+1], dtype=np.float32)
901+
902+
xx[0:nt, 0:nz] = np.transpose(h5data.xsurf[0:nz, 0:nt])
903+
xx[ nt, 0:nz] = xx[0 , 0:nz]
904+
xx[0:nt, nz] = xx[0:nt, 0 ]
905+
xx[ nt, nz] = xx[0 , 0 ]
906+
907+
yy[0:nt, 0:nz] = np.transpose(h5data.ysurf[0:nz, 0:nt])
908+
yy[ nt, 0:nz] = yy[0 , 0:nz]
909+
yy[0:nt, nz] = yy[0:nt, 0 ]
910+
yy[ nt, nz] = yy[0 , 0 ]
911+
912+
zz[0:nt, 0:nz] = np.transpose(h5data.zsurf[0:nz, 0:nt])
913+
zz[ nt, 0:nz] = zz[0 , 0:nz]
914+
zz[0:nt, nz] = zz[0:nt, 0 ]
915+
zz[ nt, nz] = zz[0 , 0 ]
916+
917+
#plot
918+
if plottype == 'surf':
919+
mlab.mesh(xx, yy, zz, color=color)
920+
elif plottype == 'targetBn':
921+
for i in range(nz+1):
922+
ator = (i+0.5)*2*np.pi/nz #zeta
923+
for j in range(nt+1):
924+
apol = (j+0.5)*2*np.pi/nt #theta
925+
tmpBn = h5data.target_Bmnc*np.cos(h5data.Bmnim*apol-h5data.Bmnin*ator) \
926+
+ h5data.target_Bmns*np.sin(h5data.Bmnim*apol-h5data.Bmnin*ator)
927+
928+
Bn[j,i] = np.sum(tmpBn)
929+
mlab.mesh(xx, yy, zz, scalars=Bn, colormap='RdBu')
930+
axes = mlab.colorbar( orientation='vertical', label_fmt='%.1e')
931+
axes.title_text_property.font_family = 'times'
932+
axes.title_text_property.font_size = 6
933+
axes.title_text_property.italic = False
934+
axes.label_text_property.font_family = 'times'
935+
axes.label_text_property.font_size = 100
936+
axes.label_text_property.italic = False
937+
elif plottype == 'realizedBn':
938+
Bn[0:nt, 0:nz] = np.transpose(h5data.Bn[0:nz, 0:nt])
939+
Bn[ nt, 0:nz] = Bn[0 , 0:nz]
940+
Bn[0:nt, nz] = Bn[0:nt, 0 ]
941+
Bn[ nt, nz] = Bn[0 , 0 ]
942+
mlab.mesh(xx, yy, zz, scalars=Bn, colormap='RdBu')
943+
axes = mlab.colorbar(orientation='vertical', label_fmt='%.1e')
944+
axes.title_text_property.font_family = 'times'
945+
axes.title_text_property.font_size = 6
946+
axes.title_text_property.italic = False
947+
axes.label_text_property.font_family = 'times'
948+
axes.label_text_property.font_size = 100
949+
axes.label_text_property.italic = False
950+
else :
951+
raise NameError("No such option!")
952+
print "plottype = surf/targetBn/realizedBn"
953+
return

0 commit comments

Comments
 (0)