Skip to content

Commit b946dcd

Browse files
pguenthpguenth
committed
Many changes; thesis submission
1 parent c99fc19 commit b946dcd

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

43 files changed

+5443
-359
lines changed

‎.gitignore

+4
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,8 @@ build/
33
lib/
44
out/
55
tests/
6+
testcache/
7+
tests_old/
68
pickle/
79
*/.undodir/
810
__Tagbar__.*
@@ -12,5 +14,7 @@ docs/doxy/xml
1214
tests_old/
1315
testcache/
1416
figures-old/
17+
figures-backup/
1518
figures/
1619
log/
20+
cysignals_crash_logs/

‎Makefile

+1-1
Original file line numberDiff line numberDiff line change
@@ -54,7 +54,7 @@ bin: $(CXX_OBJS)
5454

5555
clib: dirs libbatch.so libbatch.a
5656

57-
cython: dirs libbatch.a
57+
cython: dirs
5858
echo \$asdf
5959
python setup.py build_ext --build-lib $(LIB_DIR) --build-temp $(BUILD_DIR)
6060

‎achterberg-fig2-ln-correct.py

+146
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,146 @@
1+
import sys
2+
sys.path.insert(0, 'lib')
3+
sys.path.insert(0, 'src/evaluation')
4+
from pybatch.special.kruells import *
5+
import proplot as pplt
6+
import logging
7+
import chains
8+
import formats
9+
from node.cache import PickleNodeCache
10+
11+
from powerlawseries import *
12+
13+
pplt.rc.update({
14+
'text.usetex' : True,
15+
})
16+
17+
logging.basicConfig(level=logging.INFO,
18+
format='%(asctime)s - %(name)s - %(levelname)s - %(message)s')
19+
20+
def cb(this_param, param, other_param):
21+
dt = other_param['alpha'] / param['V']
22+
q = param['V'] * other_param['Ldiff']
23+
Ls = other_param['Ldiff'] * this_param['epsilon']
24+
batch_param = param | {'q' : q, 'dt' : dt, 'Ls' : Ls}
25+
label_fmt_fields = this_param | other_param
26+
return {'param': batch_param, 'label_fmt_fields': label_fmt_fields}
27+
28+
param = {
29+
'r' : 4,
30+
'x0' : 0,
31+
'y0' : 0,
32+
't_inj' : 1, #0.1, # 0.2
33+
'Tmax' : 3500,
34+
'V' : 1,
35+
}
36+
other_param = {
37+
'alpha' : 0.05,
38+
'Ldiff' : 1,
39+
}
40+
41+
name = 'achterberg_fig2_ln_correct'
42+
cache = PickleNodeCache('pickle', name)
43+
add_conf = [(1, 0, np.inf)]
44+
chain = chains.get_chain_powerlaw_datapoint(PyBatchAchterberg2, cache, np.inf, lambda c: c['label_fmt_fields']['epsilon'], histo_opts={'bin_count' : 40}, negate_index=True, log_bins=False, additional_confine_ranges=add_conf)
45+
#chain.map_tree(lambda n: setattr(n, 'ignore_cache', True), "values")
46+
chain.map_tree(lambda n: n.set(ln_x=True) , "pl")
47+
chain.map_tree(lambda n: print(n.plot_on), "pl")
48+
49+
#var = PowerlawSeriesVariable('\\epsilon', 'epsilon', [0.01, 0.02, 0.04, 0.05, 0.08, 0.1, 0.2, 0.3, 0.4, 0.5, 0.8, 1])
50+
#var = PowerlawSeriesVariable('\\epsilon', 'epsilon', [0.01, 0.04, 0.08, 0.2, 0.4, 0.8, 1.0]) #[0.01, 0.05, 0.1, 0.8])[float(sys.argv[1])])
51+
var = PowerlawSeriesVariable('\\epsilon', 'epsilon', [float(sys.argv[1])])
52+
#var = PowerlawSeriesVariable('\\epsilon', 'epsilon', [0.01 * int(sys.argv[1])])
53+
pls = PowerlawSeries(chain, var, cb,
54+
callback_kwargs={'param': param, 'other_param': other_param}
55+
)
56+
57+
axs_format = SliceDict()
58+
axs_format[2] = dict(xlabel='ln p/p_inj', xscale='linear')
59+
axs_format[3] = dict(xlabel='ln p/p_inj', xscale='linear')
60+
axs_format[4] = dict(xlabel='$\\epsilon$', ylabel='Powerlaw index $s$', xscale='log')
61+
axs_format[slice(None, None, None)] = dict(toplabels=('Spatial distribution', 'Momentum distribution', 'Powerlaw indizes'), leftlabels=('CES', 'KPPC'))
62+
nfmt = NodeFigureFormat(base=formats.histsandpowerlaw2, fig_legend_kw=None, axs_format=axs_format, legends_kw={0: None, 1: None, 2: None, 4:{}})
63+
nfig = NodeFigure(nfmt)
64+
nfig.format(suptitle='Fig. 2 of Achterberg/Schure')
65+
chain_x, chain_p = pls.histogram_chains
66+
67+
datarow = pls.datarow_chain
68+
datarow.set(label='Cauchy-Euler scheme')
69+
70+
71+
def mod_scheme(pls, batch_cls, name, name_human):
72+
memo = {}
73+
chain_x_old, chain_p_old = pls.histogram_chains
74+
datarow_old = pls.datarow_chain
75+
chain_x = chain_x_old.copy(name, last_kwargs={'batch_cls': batch_cls}, memo=memo)
76+
chain_p = chain_p_old.copy(name, last_kwargs={'batch_cls': batch_cls}, memo=memo)
77+
datarow = datarow_old.copy(name, last_parents={'batch_cls': batch_cls}, memo=memo)
78+
datarow.set(label=name_human)
79+
return datarow, chain_x, chain_p
80+
81+
datarow_semiimplicit, chain_x_semiimplicit, chain_p_semiimplicit = mod_scheme(pls, PyBatchAchterberg2SemiImplicit2, 'semiimplicit', 'Semi-implicit')
82+
datarow_implicit, chain_x_implicit, chain_p_implicit = mod_scheme(pls, PyBatchAchterberg2Implicit, 'implicit', 'Implicit')
83+
datarow_kppc, chain_x_kppc, chain_p_kppc = mod_scheme(pls, PyBatchAchterberg2KPPC, 'kppc', 'KPPC')
84+
85+
#memo = {}
86+
#secondorder_cls = PyBatchAchterberg2SecondOrder
87+
#chain_x_2nd = chain_x.copy("2ndorder", last_kwargs={'batch_cls': secondorder_cls}, memo=memo)
88+
#chain_p_2nd = chain_p.copy("2ndorder", last_kwargs={'batch_cls': secondorder_cls}, memo=memo)
89+
#datarow_2nd = datarow.copy("2ndorder", last_parents={'batch_cls': secondorder_cls}, memo=memo)
90+
#datarow_2nd.set(label='2nd order scheme')
91+
#
92+
#memo = {}
93+
#secondorder2_cls = PyBatchAchterberg2SecondOrder2
94+
#chain_x_2nd2 = chain_x.copy("2ndorder2", last_kwargs={'batch_cls': secondorder2_cls}, memo=memo)
95+
#chain_p_2nd2 = chain_p.copy("2ndorder2", last_kwargs={'batch_cls': secondorder2_cls}, memo=memo)
96+
#datarow_2nd2 = datarow.copy("2ndorder2", last_parents={'batch_cls': secondorder2_cls}, memo=memo)
97+
#datarow_2nd2.set(label='2nd order scheme (vec)')
98+
99+
#memo = {}
100+
#semiimplicit_cls = PyBatchAchterberg2SemiImplicit
101+
#chain_x_semiimplicit = chain_x.copy("semiimplicit", last_kwargs={'batch_cls': semiimplicit_cls}, memo=memo)
102+
#chain_p_semiimplicit = chain_p.copy("semiimplicit", last_kwargs={'batch_cls': semiimplicit_cls}, memo=memo)
103+
#datarow_semiimplicit = datarow.copy("semiimplicit", last_parents={'batch_cls': semiimplicit_cls}, memo=memo)
104+
#datarow_semiimplicit.set(label='semiimplicit scheme (vec)')
105+
106+
#nfig.add(chain_x, 0, plot_on='spectra')
107+
#nfig.add(chain_p, 2, plot_on='spectra')
108+
nfig.add(datarow, 4)
109+
110+
#nfig.add(chain_x_kppc, 1, plot_on='spectra')
111+
#nfig.add(chain_p_kppc, 3, plot_on='spectra')
112+
nfig.add(datarow_kppc, 4)
113+
#
114+
nfig.add(chain_x_implicit, 0, plot_on='spectra')
115+
nfig.add(chain_p_implicit, 2, plot_on='spectra')
116+
nfig.add(datarow_implicit, 4)
117+
118+
comparison_function = lambda eps : 1 + 0.924 * eps + 0.095 * eps**2
119+
comparison_node = FunctionNode("fun", callback=comparison_function, plot=True)
120+
nfig.add(comparison_node, 4)
121+
122+
#nfig.add(chain_x_2nd2, 1, plot_on='spectra')
123+
#nfig.add(chain_p_2nd2, 3, plot_on='spectra')
124+
#nfig.add(datarow_2nd2, 4)
125+
126+
#nfig.add(chain_x_2nd, 1, plot_on='spectra')
127+
#nfig.add(chain_p_2nd, 3, plot_on='spectra')
128+
#nfig.add(datarow_2nd, 4)
129+
130+
nfig.add(chain_x_semiimplicit, 0, plot_on='spectra')
131+
nfig.add(chain_p_semiimplicit, 2, plot_on='spectra')
132+
nfig.add(datarow_semiimplicit, 4)
133+
134+
135+
nfig.pad(0.2, 4)
136+
nfig[2].legend()
137+
nfig[3].legend()
138+
#nfig.show_nodes("achterberg_tree.pdf")
139+
nfig.savefig("figures/{}.pdf".format(name))
140+
141+
#nfig = NodeFigure(NodeFigureFormat(subplots={'ncols' : 1}, ax_format={'xscale': 'log', 'xlabel': 'streaming velocity / diffusivity', 'ylabel': 'powerlaw index'}, legend_kw={'ncols': 1}))
142+
#nfig.add(pls.datarow_chain)
143+
#nfig.add(datarow_kppc)
144+
#nfig.add(datarow_implicit)
145+
#nfig.add(comparison_node)
146+
#nfig.savefig("figures/{}-simple.pdf".format(name))

‎achterberg-fig2-ln.py

+86-14
Original file line numberDiff line numberDiff line change
@@ -28,8 +28,8 @@ def cb(this_param, param, other_param):
2828
'r' : 4,
2929
'x0' : 0,
3030
'y0' : 0,
31-
't_inj' : 0.2,
32-
'Tmax' : 10000,
31+
't_inj' : 0.1, # 0.2
32+
'Tmax' : 3500,
3333
'V' : 1,
3434
'Ls' : 1,
3535
}
@@ -39,12 +39,14 @@ def cb(this_param, param, other_param):
3939

4040
name = 'achterberg_fig2_ln'
4141
cache = PickleNodeCache('pickle', name)
42-
chain = chains.get_chain_powerlaw_datapoint(PyBatchAchterberg2, cache, np.inf, lambda c: c['label_fmt_fields']['epsilon'], histo_opts={'bin_count' : 50}, negate_index=True, log_bins=False)
43-
chain.map_tree(lambda n: setattr(n, 'ignore_cache', True), "values")
44-
chain.map_tree(lambda n: n.set(ln_x=True), "pl")
42+
add_conf = [(1, 0, np.inf)]
43+
chain = chains.get_chain_powerlaw_datapoint(PyBatchAchterberg2, cache, np.inf, lambda c: c['label_fmt_fields']['epsilon'], histo_opts={'bin_count' : 40}, negate_index=True, log_bins=False, additional_confine_ranges=add_conf)
44+
#chain.map_tree(lambda n: setattr(n, 'ignore_cache', True), "values")
45+
chain.map_tree(lambda n: n.set(ln_x=True) , "pl")
46+
chain.map_tree(lambda n: print(n.plot_on), "pl")
4547

4648
#var = PowerlawSeriesVariable('\\epsilon', 'epsilon', [0.01, 0.02, 0.04, 0.05, 0.08, 0.1, 0.2, 0.3, 0.4, 0.5, 0.8, 1])
47-
var = PowerlawSeriesVariable('\\epsilon', 'epsilon', [0.01, 0.02, 0.04, 0.08, 0.2, 0.4, 0.8, 1])
49+
var = PowerlawSeriesVariable('\\epsilon', 'epsilon', [0.01, 0.04, 0.08, 0.2, 0.4, 0.8, 1.0]) #[0.01, 0.05, 0.1, 0.8])[float(sys.argv[1])])
4850
#var = PowerlawSeriesVariable('\\epsilon', 'epsilon', [0.01 * int(sys.argv[1])])
4951
pls = PowerlawSeries(chain, var, cb,
5052
callback_kwargs={'param': param, 'other_param': other_param}
@@ -60,22 +62,92 @@ def cb(this_param, param, other_param):
6062
nfig.format(suptitle='Fig. 2 of Achterberg/Schure')
6163
chain_x, chain_p = pls.histogram_chains
6264

65+
datarow = pls.datarow_chain
66+
datarow.set(label='Cauchy-Euler scheme')
67+
6368
memo = {}
6469
kppc_cls = PyBatchAchterberg2KPPC
6570
chain_x_kppc = chain_x.copy("kppc", last_kwargs={'batch_cls': kppc_cls}, memo=memo)
6671
chain_p_kppc = chain_p.copy("kppc", last_kwargs={'batch_cls': kppc_cls}, memo=memo)
67-
datarow = pls.datarow_chain
68-
datarow.set(label='CES')
6972
datarow_kppc = datarow.copy("kppc", last_parents={'batch_cls': kppc_cls}, memo=memo)
70-
datarow_kppc.set(label='KPPC')
73+
datarow_kppc.set(label='Predictor-corrector scheme')
74+
75+
memo = {}
76+
implicit_cls = PyBatchAchterberg2Implicit
77+
chain_x_implicit = chain_x.copy("implicit", last_kwargs={'batch_cls': implicit_cls}, memo=memo)
78+
chain_p_implicit = chain_p.copy("implicit", last_kwargs={'batch_cls': implicit_cls}, memo=memo)
79+
datarow_implicit = datarow.copy("implicit", last_parents={'batch_cls': implicit_cls}, memo=memo)
80+
datarow_implicit.set(label='Implicit Euler scheme')
81+
82+
#memo = {}
83+
#secondorder_cls = PyBatchAchterberg2SecondOrder
84+
#chain_x_2nd = chain_x.copy("2ndorder", last_kwargs={'batch_cls': secondorder_cls}, memo=memo)
85+
#chain_p_2nd = chain_p.copy("2ndorder", last_kwargs={'batch_cls': secondorder_cls}, memo=memo)
86+
#datarow_2nd = datarow.copy("2ndorder", last_parents={'batch_cls': secondorder_cls}, memo=memo)
87+
#datarow_2nd.set(label='2nd order scheme')
88+
#
89+
#memo = {}
90+
#secondorder2_cls = PyBatchAchterberg2SecondOrder2
91+
#chain_x_2nd2 = chain_x.copy("2ndorder2", last_kwargs={'batch_cls': secondorder2_cls}, memo=memo)
92+
#chain_p_2nd2 = chain_p.copy("2ndorder2", last_kwargs={'batch_cls': secondorder2_cls}, memo=memo)
93+
#datarow_2nd2 = datarow.copy("2ndorder2", last_parents={'batch_cls': secondorder2_cls}, memo=memo)
94+
#datarow_2nd2.set(label='2nd order scheme (vec)')
95+
96+
#memo = {}
97+
#semiimplicit_cls = PyBatchAchterberg2SemiImplicit
98+
#chain_x_semiimplicit = chain_x.copy("semiimplicit", last_kwargs={'batch_cls': semiimplicit_cls}, memo=memo)
99+
#chain_p_semiimplicit = chain_p.copy("semiimplicit", last_kwargs={'batch_cls': semiimplicit_cls}, memo=memo)
100+
#datarow_semiimplicit = datarow.copy("semiimplicit", last_parents={'batch_cls': semiimplicit_cls}, memo=memo)
101+
#datarow_semiimplicit.set(label='semiimplicit scheme (vec)')
102+
103+
memo = {}
104+
semiimplicit_cls = PyBatchAchterberg2SemiImplicit2
105+
chain_x_semiimplicit2 = chain_x.copy("semiimplicit2", last_kwargs={'batch_cls': semiimplicit_cls}, memo=memo)
106+
chain_p_semiimplicit2 = chain_p.copy("semiimplicit2", last_kwargs={'batch_cls': semiimplicit_cls}, memo=memo)
107+
datarow_semiimplicit2 = datarow.copy("semiimplicit2", last_parents={'batch_cls': semiimplicit_cls}, memo=memo)
108+
datarow_semiimplicit2.set(label='Semiimplicit scheme')
109+
110+
#nfig.add(chain_x, 0, plot_on='spectra')
111+
#nfig.add(chain_p, 2, plot_on='spectra')
112+
nfig.add(datarow, 4)
71113

72-
nfig.add(chain_x, 0, plot_on='spectra')
73-
nfig.add(chain_p, 2, plot_on='spectra')
74-
nfig.add(chain_x_kppc, 1, plot_on='spectra')
75-
nfig.add(chain_p_kppc, 3, plot_on='spectra')
76-
nfig.add(pls.datarow_chain, 4)
114+
#nfig.add(chain_x_kppc, 1, plot_on='spectra')
115+
#nfig.add(chain_p_kppc, 3, plot_on='spectra')
77116
nfig.add(datarow_kppc, 4)
117+
#
118+
nfig.add(chain_x_implicit, 0, plot_on='spectra')
119+
nfig.add(chain_p_implicit, 2, plot_on='spectra')
120+
nfig.add(datarow_implicit, 4)
121+
122+
comparison_function = lambda eps : 1 + 0.924 * eps + 0.095 * eps**2
123+
comparison_node = FunctionNode("fun", callback=comparison_function, plot=True)
124+
nfig.add(comparison_node, 4)
125+
126+
#nfig.add(chain_x_2nd2, 1, plot_on='spectra')
127+
#nfig.add(chain_p_2nd2, 3, plot_on='spectra')
128+
#nfig.add(datarow_2nd2, 4)
129+
130+
#nfig.add(chain_x_2nd, 1, plot_on='spectra')
131+
#nfig.add(chain_p_2nd, 3, plot_on='spectra')
132+
#nfig.add(datarow_2nd, 4)
133+
134+
#nfig.add(chain_x_semiimplicit, 0, plot_on='spectra')
135+
#nfig.add(chain_p_semiimplicit, 2, plot_on='spectra')
136+
#nfig.add(datarow_semiimplicit, 4)
137+
138+
nfig.add(chain_x_semiimplicit2, 1, plot_on='spectra')
139+
nfig.add(chain_p_semiimplicit2, 3, plot_on='spectra')
140+
nfig.add(datarow_semiimplicit2, 4)
141+
78142
nfig.pad(0.2, 4)
143+
nfig[2].legend()
144+
nfig[3].legend()
79145
#nfig.show_nodes("achterberg_tree.pdf")
80146
nfig.savefig("figures/{}.pdf".format(name))
81147

148+
#nfig = NodeFigure(NodeFigureFormat(subplots={'ncols' : 1}, ax_format={'xscale': 'log', 'xlabel': 'streaming velocity / diffusivity', 'ylabel': 'powerlaw index'}, legend_kw={'ncols': 1}))
149+
#nfig.add(pls.datarow_chain)
150+
#nfig.add(datarow_kppc)
151+
#nfig.add(datarow_implicit)
152+
#nfig.add(comparison_node)
153+
#nfig.savefig("figures/{}-simple.pdf".format(name))

‎boundary_influence.py

+26-1
Original file line numberDiff line numberDiff line change
@@ -7,6 +7,7 @@
77
import chains
88
import formats
99
from node.cache import PickleNodeCache
10+
from matplotlib import ticker
1011

1112
from powerlawseries import *
1213

@@ -33,13 +34,37 @@ def cb(this_param, param):
3334
'beta_s' : 0.57,
3435
'q' : 5
3536
}
37+
# dxadv(') = 0.09975
38+
# dxdiff = 0.285
39+
# delta = 0.35
40+
# sigma = 0.811067
3641

3742

3843
name = 'boundary_11'
3944
cache = PickleNodeCache('pickle', name)
4045
chain = chains.get_chain_powerlaw_datapoint(PyBatchKruells11, cache, 2, lambda c: c['batch_param']['L'])
4146

42-
var = PowerlawSeriesVariable('L', 'L', [3, 5, 6, 8, 10, 20, 30, 50, 80, 200, 1000])
47+
var = PowerlawSeriesVariable('L', 'L', [3, 5, 6, 8, 10, 20, 30, 50, 80])#, 200, 1000])
4348
pls = PowerlawSeries(chain, var, cb, callback_kwargs={'param': param})
4449
pls.plot_datarow("figures/{}.pdf".format(name), formats.powerlaws, xlabel='$L$', xscale='log')
50+
nfig = NodeFigure(formats.powerlaws, title='Spectral steepening at free-escape boundary', xlabel='Boundary radius (from shock) $X_C$', xscale='log')
51+
nfig.add(pls.datarow_chain)
52+
nfig.pad(.2)
53+
54+
# endlevel hack for minor ticks
55+
logf = ticker.LogFormatter(labelOnlyBase=False, minor_thresholds=(10, 0.4))
56+
nfig[0].xaxis.set_minor_formatter(logf)
57+
def formatterfunc(x, pos):
58+
logfn = logf(x, pos)
59+
if logfn[0] in ['5', '7', '8', '9']:
60+
return ""
61+
else:
62+
return "$" + logf(x, pos) + "$"
63+
nfig[0].xaxis.set_minor_formatter(ticker.FuncFormatter(formatterfunc))
64+
nfig[0].xaxis.set_minor_locator(ticker.LogLocator(subs='all'))
65+
# ----------------------------
66+
67+
nfig[0].annotate('$\\Delta x_\\mathrm{diff} = 0.285,~~\\delta =0.35,~~\\sigma\\approx0.81$', (0.32, 0.18), xycoords='figure fraction', bbox=dict(boxstyle="square,pad=0.5", fc="white", ec="black", lw=1))
68+
nfig._legends_kw = {}
69+
nfig.savefig("figures/{}.pdf".format(name))
4570
pls.plot_histograms("figures/{}_histograms.pdf".format(name), formats.doublehist)

0 commit comments

Comments
 (0)