Skip to content
Snippets Groups Projects
Commit d30049a0 authored by rem's avatar rem
Browse files

restructure modules, small update tbmdata_aug 

git-svn-id: https://solps-mdsplus.aug.ipp.mpg.de/repos/TORBEAM/branches/python@532 438fab6c-3626-0410-97ee-e5509541aacd
parent af630788
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
augdemo.py
+ 19
− 19
View file @ d30049a0
......@@ -2,7 +2,7 @@
from util import *
from ctbm import *
from getECangles import *
from tbmdata_aug import getECRHanglesAUG
import matplotlib.pyplot as plt
system = 0
......@@ -13,31 +13,31 @@ power = 0.6806
# Build data from shotfile:
# import tbmdata
# dat = tbmdata.tbmd(shot)
# dat.readAll()
# m, n, eqdata = dat.eqdata(time=time)
# k, l, prdata = dat.prdata(time=time)
# intinbeam, dblinbeam = dat.default_inbeams()
# dblinbeam = dat.beamdata(time=time, gyr=systemi, beamdbl=dblinbeam)
import tbmdata_aug as tbmdata
dat = tbmdata.tbmd(shot)
dat.readAll()
m, n, eqdata = dat.eqdata(time=time)
k, l, prdata = dat.prdata(time=time)
intinbeam, dblinbeam = dat.default_inbeams()
dblinbeam = dat.beamdata(time=time, gyr=system, beamdbl=dblinbeam)
# Build eqdata from input files
m, n, eqdata = readTopfile3('augdata/topfile');
k, l, prdata = readPrdata('augdata/ne.dat', 'augdata/Te.dat');
intinbeam, dblinbeam = default_inbeams();
#m, n, eqdata = readTopfile3('augdata/topfile');
#k, l, prdata = readPrdata('augdata/ne.dat', 'augdata/Te.dat');
#intinbeam, dblinbeam = default_inbeams();
# Set beam power
dblinbeam[23] = power;
#dblinbeam[23] = power;
# Set launcher position
dblinbeam[3] = 236.1;
dblinbeam[5] = 32.;
#dblinbeam[3] = 236.1;
#dblinbeam[5] = 32.;
# AUG specifics
if system == 7 or system == 8:
dblinbeam[5] = -32.
#if system == 7 or system == 8:
# dblinbeam[5] = -32.
# read angles from shotfile
if system > 0:
pol, tor = getECRHanglesAUG(shot, system, time=time)
dblinbeam[2], dblinbeam[1] = ecrh2tbm(pol, tor, system, shot);
#if system > 0:
# pol, tor = getECRHanglesAUG(shot, system, time=time)
# dblinbeam[2], dblinbeam[1] = ecrh2tbm(pol, tor, system, shot);
# set working absorption routine (Farina)
intinbeam[11] = 1
......
ctbm.py
+ 2
− 1
View file @ d30049a0
......@@ -7,7 +7,8 @@ import re
import sys
from ctypes import *
libtbm = '/afs/ipp/u/rem/soft/libTORBEAM/lib-OUT/libtorbeam.so'
libtbm = '../lib-OUT/libtorbeam.so'
#libtbm = './libtorbeam.so'
tbm = cdll.LoadLibrary(libtbm)
def readTopfile(fname):
......
tbmdata.py → tbmdata_aug.py
+ 28
− 2
View file @ d30049a0
......@@ -18,8 +18,33 @@ try:
system = os.environ['OS'];
except KeyError:
system = 'amd64_sles11';
libecrh = '/afs/ipp/home/e/ecrh/'+system+'/lib/libaug_ecrh_setmirrors.so';
try:
libecrh = '/afs/ipp/home/e/ecrh/amd64_sles11/lib/libaug_ecrh_setmirrors.so';
ecrh = cdll.LoadLibrary(libecrh)
except:
print('Cannot find AUG ECRH library, will not run AUG live cases');
def getECRHanglesAUG(shot, jgy, time=0.0):
if jgy < 5:
sys = 'P_sy1_g'+str(jgy);
else:
sys = 'P_sy2_g'+str(jgy - 4);
sf = dd.shotfile('ECS', shot);
pol = sf.getParameter(sys, 'GPolPos').data;
tor = sf.getParameter(sys, 'GTorPos').data;
sf.close();
if jgy > 4 and time > 0:
sf = dd.shotfile('ECN', shot);
sig = sf("G"+str(jgy-4)+"POL");
if (time < sig.time[0] or time > sig.time[sig.time.size-1]):
raise Exception("Time not accessible in mirror position signal");
pol = sig.data[sig.time > time][0]/100; # [V] (=cm) -> [m]
sf.close();
return pol, tor;
def ecrh2tbm(pol_ecs, tor_ecs, jgy, nshot):
"""Converts ECRH machine parameters to TORBEAM compatible angles
......@@ -73,6 +98,7 @@ class tbmd:
self.eqdok = True
self.kk._read_pfm()
self.kk._read_scalars()
self.kk._read_grid()
else:
print('Cannot continue without equilibrium data.')
......
getECangles.py → tbmdata_tcv.py
+ 51
− 71
View file @ d30049a0
#!/usr/bin/env python
try:
import dd
except:
print('Module dd not found, cannot read AUG data from shotfiles');
import re
import sys
import numpy as np
from ctypes import *
from subprocess import Popen, PIPE
try:
libecrh = '/afs/ipp/home/e/ecrh/amd64_sles11/lib/libaug_ecrh_setmirrors.so';
ecrh = cdll.LoadLibrary(libecrh)
except:
print('Cannot find AUG ECRH library, will not run AUG live cases');
def getECRHanglesAUG(shot, jgy, time=0.0):
if jgy < 5:
sys = 'P_sy1_g'+str(jgy);
else:
sys = 'P_sy2_g'+str(jgy - 4);
sf = dd.shotfile('ECS', shot);
pol = sf.getParameter(sys, 'GPolPos').data;
tor = sf.getParameter(sys, 'GTorPos').data;
sf.close();
if jgy > 4 and time > 0:
sf = dd.shotfile('ECN', shot);
sig = sf("G"+str(jgy-4)+"POL");
if (time < sig.time[0] or time > sig.time[sig.time.size-1]):
raise Exception("Time not accessible in mirror position signal");
pol = sig.data[sig.time > time][0]/100; # [V] (=cm) -> [m]
sf.close();
return pol, tor;
def ecrh2tbm(pol_ecs, tor_ecs, jgy, nshot):
"""Converts ECRH machine parameters to TORBEAM compatible angles
Input can be directly from values read from shotfile.
System ECRH1 has parameters GPolPos and GTorPos in parameter set
P_sy1_g which don't change during a shot.
System ECRH2 has parameters beta and time-varying poloidal mirror
which can be read from ECN.
"""
if nshot > 27400:
datum = 0;
else:
datum = 20000101;
c_err = c_int32(0);
c_sysunt = c_int32(100 + 100*(jgy / 5) + (jgy % 5) + (jgy / 5));
if jgy > 4:
c_theta = c_double(1.e3*pol_ecs); # [m] -> [mm]
else:
c_theta = c_double(pol_ecs);
c_phi = c_double(tor_ecs);
_err = byref(c_err);
_sysunt = byref(c_sysunt);
_theta = byref(c_theta);
_phi = byref(c_phi);
c_datum = c_double(datum);
_datum = byref(c_datum);
s = ecrh.setval2tp_(_err, _sysunt, _theta, _phi, _datum);
if c_err.value < 0 or c_err.value > 100:
raise Exception("Parameter problem: "+str(c_err.value));
if c_err.value > 0:
raise Exception("Calculation problem: "+str(c_err.value));
return -c_theta.value, -c_phi.value;
def cmdline(command):
process = Popen(
args=command,
......@@ -182,3 +112,53 @@ def write_inbeam_tcv(shot, time, launcher):
f.write(' xqedg = 4.000000 ! safety factor on the edge\n');
f.write('/\n');
def file2array(fname, tplidx=0):
f = open(fname);
strarr = np.array(f.readlines());
f.close();
strarr = strarr[tplidx:];
nnum = len( re.split('\s+', strarr[0].strip()) );
mnum = strarr.size;
res = np.zeros([nnum, mnum]);
res[:,:] = np.NAN;
ic = 0;
for s in strarr:
jc = 0
for i in re.split('\s+', s.strip()):
try:
res[jc][ic] = float(i);
jc += 1;
except:
jc += 1;
ic += 1;
return res
def btvac_wrap_tcv(shot,time):
script = "function [btvac] = auto_btvac(shot, time)\n" + \
"% prints to stdout value of Btor (@ 0.88 m) at time for shot\n" + \
"\n" + \
"fprintf('BEGIN OUTPUT\\n');\n" + \
"try\n" + \
" a = gdat(shot, 'b0');\n" + \
" btvac = mean(a.data(find(a.t(find(a.t < time+0.05)) > time-0.05)));\n" + \
" fprintf('>>> %f\\n', btvac);\n" + \
"catch\n" + \
" fprintf('! error from gdat, no result \\n');\n" + \
"end\n" + \
" fprintf('END OUTPUT\\n');\n" + \
"\n" + \
"end\n"
with open('auto_btvac.m','w') as f:
f.write(script);
cmd = 'matlab -nodesktop -nodisplay -nosplash -r "auto_btvac('+ str(shot)+','+str(time)+');quit;"';
print 'execute: '+cmd;
res2 = cmdline(cmd);
resL = res2[res2.index('BEGIN OUTPUT')+len('BEGIN OUTPUT\n'):res2.index('END OUTPUT')];
strarr = np.array(resL.split('\n'));
s = strarr[0][strarr[0].index('>>> ')+4:];
btvac = float(s);
cmdline('rm -f auto_btvac.m');
return btvac
tcvdemo.py
+ 3
− 3
View file @ d30049a0
......@@ -8,11 +8,11 @@ import os
import sys
import numpy as np
import scipy as sp
from itertools import chain
import eqtools as eqt
from ctbm import *
from getECangles import *
from itertools import chain
from util import btvac_wrap_tcv
from tbmdata_tcv import *
if len(sys.argv) < 2:
print('Need a TORAY dataset, identify with shot+time');
......
util.py deleted 100644 → 0
+ 0
− 65
View file @ af630788
#!/usr/bin/env python
import re
from subprocess import Popen, PIPE
import sys
import numpy as np
def cmdline(command):
process = Popen(
args=command,
stdout=PIPE,
shell=True
)
return process.communicate()[0]
def file2array(fname, tplidx=0):
f = open(fname);
strarr = np.array(f.readlines());
f.close();
strarr = strarr[tplidx:];
nnum = len( re.split('\s+', strarr[0].strip()) );
mnum = strarr.size;
res = np.zeros([nnum, mnum]);
res[:,:] = np.NAN;
ic = 0;
for s in strarr:
jc = 0
for i in re.split('\s+', s.strip()):
try:
res[jc][ic] = float(i);
jc += 1;
except:
jc += 1;
ic += 1;
return res
def btvac_wrap_tcv(shot,time):
script = "function [btvac] = auto_btvac(shot, time)\n" + \
"% prints to stdout value of Btor (@ 0.88 m) at time for shot\n" + \
"\n" + \
"fprintf('BEGIN OUTPUT\\n');\n" + \
"try\n" + \
" a = gdat(shot, 'b0');\n" + \
" btvac = mean(a.data(find(a.t(find(a.t < time+0.05)) > time-0.05)));\n" + \
" fprintf('>>> %f\\n', btvac);\n" + \
"catch\n" + \
" fprintf('! error from gdat, no result \\n');\n" + \
"end\n" + \
" fprintf('END OUTPUT\\n');\n" + \
"\n" + \
"end\n"
with open('auto_btvac.m','w') as f:
f.write(script);
cmd = 'matlab -nodesktop -nodisplay -nosplash -r "auto_btvac('+ str(shot)+','+str(time)+');quit;"';
print 'execute: '+cmd;
res2 = cmdline(cmd);
resL = res2[res2.index('BEGIN OUTPUT')+len('BEGIN OUTPUT\n'):res2.index('END OUTPUT')];
strarr = np.array(resL.split('\n'));
s = strarr[0][strarr[0].index('>>> ')+4:];
btvac = float(s);
cmdline('rm -f auto_btvac.m');
return btvac
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment