-
Notifications
You must be signed in to change notification settings - Fork 0
/
generate_obsidian_summary.py
118 lines (100 loc) · 4.04 KB
/
generate_obsidian_summary.py
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
import booz_xform as bx
from simsopt.mhd.vmec import Vmec
import argparse
import os
from simsopt._core.optimizable import load, save
from simsopt.geo.surfacerzfourier import SurfaceRZFourier
from simsopt.field import BiotSavart
import numpy as np
import matplotlib.pyplot as plt
from simsopt.field import SurfaceClassifier, \
particles_to_vtk, compute_fieldlines, LevelsetStoppingCriterion, plot_poincare_data
from simsopt.geo import CurveLength
import pickle
# Read command line arguments
parser = argparse.ArgumentParser()
parser.add_argument("--opt", dest="opt", default=None)
parser.add_argument("--MN", dest="MN", default=None)
# Prepare args
args = parser.parse_args()
pp = f'runs/opt_{args.opt}/M=N={args.MN}'
name = f'OPT_{args.opt}_M=N={args.MN}'
this_path = os.path.join(os.getcwd(), pp)
figure_path = os.path.join( this_path, 'figure')
obsidian_path = '/Users/antoinebaillod/Library/Mobile Documents/iCloud~md~obsidian/Documents/Research/CORE/RUNS'
os.makedirs(figure_path, exist_ok=True)
os.chdir(this_path)
v = Vmec('input.final')
bs = load(os.path.join(this_path, 'coils/bs_output.json'))
v.run()
# Run boozer
b = bx.Booz_xform()
b.read_wout(v.output_file)
b.mboz = 48
b.nboz = 48
b.run()
b.write_boozmn("boozmn.nc")
s = b.s_b
# Evaluate figures of merit
bmnc = b.bmnc_b
xm = b.xm_b
xn = b.xn_b
f = np.sqrt(np.sum(bmnc[xn!=0,:]**2,axis=0)/np.sum(bmnc**2,axis=0))
fqs = np.mean(f)
mean_iota = np.mean(b.iota)
volume = v.boundary.volume()
major_radius = v.boundary.major_radius()
minor_radius = v.boundary.minor_radius()
aspect_ratio = v.boundary.aspect_ratio()
il_current = np.abs(bs.coils[0].current.get_value())
pf_current = np.abs(bs.coils[2].current.get_value())
ll = CurveLength( bs.coils[0].curve )
il_length = ll.J()
fname = name + '.md'
with open( os.path.join(obsidian_path, fname), 'w') as file:
file.write("---\n")
file.write("tags:\n")
file.write(" - csx\n")
file.write(" - vmec_combined_approach\n")
file.write(" - optimization/run\n")
file.write("aliases:\n")
file.write(f" - OPT_{args.opt}\n")
file.write(f'name: "{name}"\n')
file.write(f'mean iota: "{mean_iota:.2E}"\n')
file.write(f'fQS: "{fqs:.2E}"\n')
file.write(f'volume: "{volume:.2E}"\n')
file.write(f'rmaj: "{major_radius:.2E}"\n')
file.write(f'rmin: "{minor_radius:.2E}"\n')
file.write(f'aspect: "{aspect_ratio:.2E}"\n')
file.write(f'IL current: "{il_current:.2E}"\n')
file.write(f'PF current: "{pf_current:.2E}"\n')
file.write(f'IL length: "{il_length:.2E}"\n')
file.write("---\n\n\n")
file.write("## Notes\n\n\n")
file.write("## Inputs\n\n")
file.write("```\n")
with open( os.path.join(this_path, 'input.txt'), 'r' ) as input_file:
for line in input_file:
file.write(line)
file.write("```\n\n\n")
file.write("## Figures\n")
file.write("### 3D plot of coils and surface\n")
file.write("### Poincare section\n")
file.write(f"})\n")
file.write("### |B| on plasma boundary\n")
file.write(f"})\n")
file.write("### $\mathbf{B}\cdot\hat{n}$ error (normalized)\n")
file.write(f"})\n")
file.write("### Profiles\n")
file.write(f"})\n")
file.write(f"})\n")
# with open( os.path.join(this_path, 'figures_of_merit.txt'), 'w' ) as file:
# file.write(f'Mean iota = {mean_iota:.2E}\n')
# file.write(f'f_quasisymmetry = {fqs:.2E}\n')
# file.write(f'Volume = {volume:.2E}\n')
# file.write(f'Major radius = {major_radius:.2E}\n')
# file.write(f'Minor radius = {minor_radius:.2E}\n')
# file.write(f'Aspect ratio = {aspect_ratio:.2E}\n')
# file.write(f'IL coil current = {il_current:.2E}\n')
# file.write(f'PF coil current = {pf_current:.2E}\n')
# file.write(f'IL coil length = {il_length:.2E}\n')