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
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
1284
1285
1286
1287
1288
1289
1290
1291
1292
1293
1294
1295
1296
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
1309
1310
1311
1312
|
<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN"
"http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
<html xmlns="http://www.w3.org/1999/xhtml">
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
<title>Worked Examples — XMDS2 2.2.2 documentation</title>
<link rel="stylesheet" href="_static/default.css" type="text/css" />
<link rel="stylesheet" href="_static/pygments.css" type="text/css" />
<script type="text/javascript">
var DOCUMENTATION_OPTIONS = {
URL_ROOT: './',
VERSION: '2.2.2',
COLLAPSE_INDEX: false,
FILE_SUFFIX: '.html',
HAS_SOURCE: true
};
</script>
<script type="text/javascript" src="_static/jquery.js"></script>
<script type="text/javascript" src="_static/underscore.js"></script>
<script type="text/javascript" src="_static/doctools.js"></script>
<script type="text/javascript" src="http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS_HTML,http://www.xmds.org/_static/mathjax-use-tex-fonts.js"></script>
<link rel="shortcut icon" href="_static/xmds_favicon.ico"/>
<link rel="top" title="XMDS2 2.2.2 documentation" href="index.html" />
<link rel="next" title="Reference section" href="reference_index.html" />
<link rel="prev" title="Quickstart Tutorial" href="tutorial.html" />
</head>
<body>
<div class="related">
<h3>Navigation</h3>
<ul>
<li class="right" style="margin-right: 10px">
<a href="genindex.html" title="General Index"
accesskey="I">index</a></li>
<li class="right" >
<a href="reference_index.html" title="Reference section"
accesskey="N">next</a> |</li>
<li class="right" >
<a href="tutorial.html" title="Quickstart Tutorial"
accesskey="P">previous</a> |</li>
<li><a href="index.html">XMDS2 2.2.2 documentation</a> »</li>
</ul>
</div>
<div class="document">
<div class="documentwrapper">
<div class="bodywrapper">
<div class="body">
<div class="section" id="worked-examples">
<span id="workedexamples"></span><span id="index-0"></span><h1>Worked Examples<a class="headerlink" href="#worked-examples" title="Permalink to this headline">¶</a></h1>
<p>One of the best ways to learn XMDS2 is to see several illustrative examples. Here are a set of example scripts and explanations of the code, which will be a good way to get started. As an instructional aid, they are meant to be read sequentially, but the adventurous could try starting with one that looked like a simulation they wanted to run, and adapt for their own purposes.</p>
<blockquote>
<div><p><a class="reference internal" href="#nonlinearschrodingerequation"><em>The nonlinear Schrödinger equation</em></a> (partial differential equation)</p>
<p><a class="reference internal" href="#kubo"><em>Kubo Oscillator</em></a> (stochastic differential equations)</p>
<p><a class="reference internal" href="#fibre"><em>Fibre Noise</em></a> (stochastic partial differential equation using parallel processing)</p>
<p><a class="reference internal" href="#integerdimensionexample"><em>Integer Dimensions</em></a> (integer dimensions)</p>
<p><a class="reference internal" href="#wignerarguments"><em>Wigner Function</em></a> (two dimensional PDE using parallel processing, passing arguments in at run time)</p>
<p><a class="reference internal" href="#groundstatebec"><em>Finding the Ground State of a BEC (continuous renormalisation)</em></a> (PDE with continual renormalisation - computed vectors, filters, breakpoints)</p>
<p><a class="reference internal" href="#hermitegaussgroundstatebec"><em>Finding the Ground State of a BEC again</em></a> (Hermite-Gaussian basis)</p>
<p><a class="reference internal" href="#dmultistatese"><em>Multi-component Schrödinger equation</em></a> (combined integer and continuous dimensions with matrix multiplication, aliases)</p>
</div></blockquote>
<p>All of these scripts are available in the included “examples” folder, along with more examples that demonstrate other tricks. Together, they provide starting points for a huge range of different simulations.</p>
<div class="section" id="the-nonlinear-schrodinger-equation">
<span id="nonlinearschrodingerequation"></span><span id="index-1"></span><h2>The nonlinear Schrödinger equation<a class="headerlink" href="#the-nonlinear-schrodinger-equation" title="Permalink to this headline">¶</a></h2>
<p>This worked example will show a range of new features that can be used in an <strong>XMDS2</strong> script, and we will also examine our first partial differential equation. We will take the one dimensional nonlinear Schrödinger equation, which is a common nonlinear wave equation. The equation describing this problem is:</p>
<div class="math">
\[\frac{\partial \phi}{\partial \xi} = \frac{i}{2}\frac{\partial^2 \phi}{\partial \tau^2} - \Gamma(\tau)\phi+i|\phi|^2 \phi\]</div>
<p>where <span class="math">\(\phi\)</span> is a complex-valued field, and <span class="math">\(\Gamma(\tau)\)</span> is a <span class="math">\(\tau\)</span>-dependent damping term. Let us look at an XMDS2 script that integrates this equation, and then examine it in detail.</p>
<div class="highlight-xpdeint"><div class="highlight"><pre><span class="nt"><simulation</span> <span class="na">xmds-version=</span><span class="s">"2"</span><span class="nt">></span>
<span class="nt"><name></span>nlse<span class="nt"></name></span>
<span class="nt"><author></span>Joe Hope<span class="nt"></author></span>
<span class="nt"><description></span>
The nonlinear Schrodinger equation in one dimension,
which is a simple partial differential equation.
We introduce several new features in this script.
<span class="nt"></description></span>
<span class="nt"><features></span>
<span class="nt"><benchmark</span> <span class="nt">/></span>
<span class="nt"><bing</span> <span class="nt">/></span>
<span class="nt"><fftw</span> <span class="na">plan=</span><span class="s">"patient"</span> <span class="nt">/></span>
<span class="nt"><openmp</span> <span class="nt">/></span>
<span class="nt"><auto_vectorise</span> <span class="nt">/></span>
<span class="nt"><globals></span>
<span class="cp"><![CDATA[</span>
<span class="k">const</span> <span class="kt">double</span> <span class="n">energy</span> <span class="o">=</span> <span class="mi">4</span><span class="p">;</span>
<span class="k">const</span> <span class="kt">double</span> <span class="n">vel</span> <span class="o">=</span> <span class="mf">0.3</span><span class="p">;</span>
<span class="k">const</span> <span class="kt">double</span> <span class="n">hwhm</span> <span class="o">=</span> <span class="mf">1.0</span><span class="p">;</span>
<span class="cp">]]></span>
<span class="nt"></globals></span>
<span class="nt"></features></span>
<span class="nt"><geometry></span>
<span class="nt"><propagation_dimension></span> xi <span class="nt"></propagation_dimension></span>
<span class="nt"><transverse_dimensions></span>
<span class="nt"><dimension</span> <span class="na">name=</span><span class="s">"tau"</span> <span class="na">lattice=</span><span class="s">"128"</span> <span class="na">domain=</span><span class="s">"(-6, 6)"</span> <span class="nt">/></span>
<span class="nt"></transverse_dimensions></span>
<span class="nt"></geometry></span>
<span class="nt"><vector</span> <span class="na">name=</span><span class="s">"wavefunction"</span> <span class="na">type=</span><span class="s">"complex"</span> <span class="na">dimensions=</span><span class="s">"tau"</span><span class="nt">></span>
<span class="nt"><components></span> phi <span class="nt"></components></span>
<span class="nt"><initialisation></span>
<span class="cp"><![CDATA[</span>
<span class="k">const</span> <span class="kt">double</span> <span class="n">w0</span> <span class="o">=</span> <span class="n">hwhm</span><span class="o">*</span><span class="n">sqrt</span><span class="p">(</span><span class="mi">2</span><span class="o">/</span><span class="n">log</span><span class="p">(</span><span class="mi">2</span><span class="p">));</span>
<span class="k">const</span> <span class="kt">double</span> <span class="n">amp</span> <span class="o">=</span> <span class="n">sqrt</span><span class="p">(</span><span class="n">energy</span><span class="o">/</span><span class="n">w0</span><span class="o">/</span><span class="n">sqrt</span><span class="p">(</span><span class="n">M_PI</span><span class="o">/</span><span class="mi">2</span><span class="p">));</span>
<span class="n">phi</span> <span class="o">=</span> <span class="n">amp</span><span class="o">*</span><span class="n">exp</span><span class="p">(</span><span class="o">-</span><span class="n">tau</span><span class="o">*</span><span class="n">tau</span><span class="o">/</span><span class="n">w0</span><span class="o">/</span><span class="n">w0</span><span class="p">)</span><span class="o">*</span><span class="n">exp</span><span class="p">(</span><span class="kc">i</span><span class="o">*</span><span class="n">vel</span><span class="o">*</span><span class="n">tau</span><span class="p">);</span>
<span class="cp">]]></span>
<span class="nt"></initialisation></span>
<span class="nt"></vector></span>
<span class="nt"><vector</span> <span class="na">name=</span><span class="s">"dampingVector"</span> <span class="na">type=</span><span class="s">"real"</span><span class="nt">></span>
<span class="nt"><components></span> Gamma <span class="nt"></components></span>
<span class="nt"><initialisation></span>
<span class="cp"><![CDATA[</span>
<span class="n">Gamma</span><span class="o">=</span><span class="mf">1.0</span><span class="o">*</span><span class="p">(</span><span class="mi">1</span><span class="o">-</span><span class="n">exp</span><span class="p">(</span><span class="o">-</span><span class="n">pow</span><span class="p">(</span><span class="n">tau</span><span class="o">*</span><span class="n">tau</span><span class="o">/</span><span class="mf">4.0</span><span class="o">/</span><span class="mf">4.0</span><span class="p">,</span><span class="mi">10</span><span class="p">)));</span>
<span class="cp">]]></span>
<span class="nt"></initialisation></span>
<span class="nt"></vector></span>
<span class="nt"><sequence></span>
<span class="nt"><integrate</span> <span class="na">algorithm=</span><span class="s">"ARK45"</span> <span class="na">interval=</span><span class="s">"20.0"</span> <span class="na">tolerance=</span><span class="s">"1e-7"</span><span class="nt">></span>
<span class="nt"><samples></span>10 100 10<span class="nt"></samples></span>
<span class="nt"><operators></span>
<span class="nt"><integration_vectors></span>wavefunction<span class="nt"></integration_vectors></span>
<span class="nt"><operator</span> <span class="na">kind=</span><span class="s">"ex"</span><span class="nt">></span>
<span class="nt"><operator_names></span>Ltt<span class="nt"></operator_names></span>
<span class="cp"><![CDATA[</span>
<span class="n">Ltt</span> <span class="o">=</span> <span class="o">-</span><span class="kc">i</span><span class="o">*</span><span class="n">ktau</span><span class="o">*</span><span class="n">ktau</span><span class="o">*</span><span class="mf">0.5</span><span class="p">;</span>
<span class="cp">]]></span>
<span class="nt"></operator></span>
<span class="cp"><![CDATA[</span>
<span class="n">dphi_dxi</span> <span class="o">=</span> <span class="n">Ltt</span><span class="p">[</span><span class="n">phi</span><span class="p">]</span> <span class="o">-</span> <span class="n">phi</span><span class="o">*</span><span class="n">Gamma</span> <span class="o">+</span> <span class="kc">i</span><span class="o">*</span><span class="nf">mod2</span><span class="p">(</span><span class="n">phi</span><span class="p">)</span><span class="o">*</span><span class="n">phi</span><span class="p">;</span>
<span class="cp">]]></span>
<span class="nt"><dependencies></span>dampingVector<span class="nt"></dependencies></span>
<span class="nt"></operators></span>
<span class="nt"></integrate></span>
<span class="nt"></sequence></span>
<span class="nt"><output></span>
<span class="nt"><sampling_group</span> <span class="na">basis=</span><span class="s">"tau"</span> <span class="na">initial_sample=</span><span class="s">"yes"</span><span class="nt">></span>
<span class="nt"><moments></span>density<span class="nt"></moments></span>
<span class="nt"><dependencies></span>wavefunction<span class="nt"></dependencies></span>
<span class="cp"><![CDATA[</span>
<span class="n">density</span> <span class="o">=</span> <span class="nf">mod2</span><span class="p">(</span><span class="n">phi</span><span class="p">);</span>
<span class="cp">]]></span>
<span class="nt"></sampling_group></span>
<span class="nt"><sampling_group</span> <span class="na">basis=</span><span class="s">"tau(0)"</span> <span class="na">initial_sample=</span><span class="s">"yes"</span><span class="nt">></span>
<span class="nt"><moments></span>normalisation<span class="nt"></moments></span>
<span class="nt"><dependencies></span>wavefunction<span class="nt"></dependencies></span>
<span class="cp"><![CDATA[</span>
<span class="n">normalisation</span> <span class="o">=</span> <span class="nf">mod2</span><span class="p">(</span><span class="n">phi</span><span class="p">);</span>
<span class="cp">]]></span>
<span class="nt"></sampling_group></span>
<span class="nt"><sampling_group</span> <span class="na">basis=</span><span class="s">"ktau(32)"</span> <span class="na">initial_sample=</span><span class="s">"yes"</span><span class="nt">></span>
<span class="nt"><moments></span>densityK<span class="nt"></moments></span>
<span class="nt"><dependencies></span>wavefunction<span class="nt"></dependencies></span>
<span class="cp"><![CDATA[</span>
<span class="n">densityK</span> <span class="o">=</span> <span class="nf">mod2</span><span class="p">(</span><span class="n">phi</span><span class="p">);</span>
<span class="cp">]]></span>
<span class="nt"></sampling_group></span>
<span class="nt"></output></span>
<span class="nt"></simulation></span>
</pre></div>
</div>
<p>Let us examine the new items in the <tt class="docutils literal"><span class="pre"><features></span></tt> element that we have demonstrated here. The existence of the <tt class="docutils literal"><span class="pre"><benchmark></span></tt> element causes the simulation to be timed. The <tt class="docutils literal"><span class="pre"><bing></span></tt> element causes the computer to make a sound upon the conclusion of the simulation. The <tt class="docutils literal"><span class="pre"><fftw></span></tt> element is used to pass options to the <a class="reference external" href="http://fftw.org">FFTW libraries for fast Fourier transforms</a>, which are needed to do spectral derivatives for the partial differential equation. Here we used the option <cite>plan=”patient”</cite>, which makes the simulation test carefully to find the fastest method for doing the FFTs. More information on possible choices can be found in the <a class="reference external" href="http://fftw.org">FFTW documentation</a>.</p>
<p>Finally, we use two tags to make the simulation run faster. The <tt class="docutils literal"><span class="pre"><auto_vectorise></span></tt> element switches on several loop optimisations that exist in later versions of the GCC compiler. The <tt class="docutils literal"><span class="pre"><openmp></span></tt> element turns on threaded parallel processing using the OpenMP standard where possible. These options are not activated by default as they only exist on certain compilers. If your code compiles with them on, then they are recommended.</p>
<p>Let us examine the <tt class="docutils literal"><span class="pre"><geometry></span></tt> element.</p>
<div class="highlight-xpdeint"><div class="highlight"><pre><span class="nt"><geometry></span>
<span class="nt"><propagation_dimension></span> xi <span class="nt"></propagation_dimension></span>
<span class="nt"><transverse_dimensions></span>
<span class="nt"><dimension</span> <span class="na">name=</span><span class="s">"tau"</span> <span class="na">lattice=</span><span class="s">"128"</span> <span class="na">domain=</span><span class="s">"(-6, 6)"</span> <span class="nt">/></span>
<span class="nt"></transverse_dimensions></span>
<span class="nt"></geometry></span>
</pre></div>
</div>
<p>This is the first example that includes a transverse dimension. We have only one dimension, and we have labelled it “tau”. It is a continuous dimension, but only defined on a grid containing 128 points (defined with the lattice variable), and on a domain from -6 to 6. The default is that transforms in continuous dimensions are fast Fourier transforms, which means that this dimension is effectively defined on a loop, and the “tau=-6” and “tau=6” positions are in fact the same. Other transforms are possible, as are discrete dimensions such as an integer-valued index, but we will leave these advanced possibilities to later examples.</p>
<p>Two vector elements have been defined in this simulation. One defines the complex-valued wavefunction “phi” that we wish to evolve. We define the transverse dimensions over which this vector is defined by the <tt class="docutils literal"><span class="pre">dimensions</span></tt> tag in the description. By default, it is defined over all of the transverse dimensions in the <tt class="docutils literal"><span class="pre"><geometry></span></tt> element, so even though we have omitted this tag for the second vector, it also assumes that the vector is defined over all of tau.</p>
<p>The second vector element contains the component “Gamma” which is a function of the transverse variable tau, as specified in the equation of motion for the field. This second vector could have been avoided in two ways. First, the function could have been written explicitly in the integrate block where it is required, but calculating it once and then recalling it from memory is far more efficient. Second, it could have been included in the “wavefunction” vector as another component, but then it would have been unnecessarily complex-valued, it would have needed an explicit derivative in the equations of motion (presumably <tt class="docutils literal"><span class="pre">dGamma_dxi</span> <span class="pre">=</span> <span class="pre">0;</span></tt>), and it would have been Fourier transformed whenever the phi component was transformed. So separating it as its own vector is far more efficient.</p>
<p>The <tt class="docutils literal"><span class="pre"><integrate></span></tt> element for a partial differential equation has some new features:</p>
<div class="highlight-xpdeint"><div class="highlight"><pre><span class="nt"><integrate</span> <span class="na">algorithm=</span><span class="s">"ARK45"</span> <span class="na">interval=</span><span class="s">"20.0"</span> <span class="na">tolerance=</span><span class="s">"1e-7"</span><span class="nt">></span>
<span class="nt"><samples></span>10 100 10<span class="nt"></samples></span>
<span class="nt"><operators></span>
<span class="nt"><integration_vectors></span>wavefunction<span class="nt"></integration_vectors></span>
<span class="nt"><operator</span> <span class="na">kind=</span><span class="s">"ex"</span><span class="nt">></span>
<span class="nt"><operator_names></span>Ltt<span class="nt"></operator_names></span>
<span class="cp"><![CDATA[</span>
<span class="n">Ltt</span> <span class="o">=</span> <span class="o">-</span><span class="kc">i</span><span class="o">*</span><span class="n">ktau</span><span class="o">*</span><span class="n">ktau</span><span class="o">*</span><span class="mf">0.5</span><span class="p">;</span>
<span class="cp">]]></span>
<span class="nt"></operator></span>
<span class="cp"><![CDATA[</span>
<span class="n">dphi_dxi</span> <span class="o">=</span> <span class="n">Ltt</span><span class="p">[</span><span class="n">phi</span><span class="p">]</span> <span class="o">-</span> <span class="n">phi</span><span class="o">*</span><span class="n">Gamma</span> <span class="o">+</span> <span class="kc">i</span><span class="o">*</span><span class="nf">mod2</span><span class="p">(</span><span class="n">phi</span><span class="p">)</span><span class="o">*</span><span class="n">phi</span><span class="p">;</span>
<span class="cp">]]></span>
<span class="nt"><dependencies></span>dampingVector<span class="nt"></dependencies></span>
<span class="nt"></operators></span>
<span class="nt"></integrate></span>
</pre></div>
</div>
<p>There are some trivial changes from the tutorial script, such as the fact that we are using the ARK45 algorithm rather than ARK89. Higher order algorithms are often better, but not always. Also, since this script has multiple output groups, we have to specify how many times each of these output groups are sampled in the <tt class="docutils literal"><span class="pre"><samples></span></tt> element, so there are three numbers there. Besides the vectors that are to be integrated, we also specify that we want to use the vector “dampingVector” during this integration. This is achieved by including the <tt class="docutils literal"><span class="pre"><dependencies></span></tt> element inside the <tt class="docutils literal"><span class="pre"><operators></span></tt> element.</p>
<p>The equation of motion as written in the CDATA block looks almost identical to our desired equation of motion, except for the term based on the second derivative, which introduces an important new concept. Inside the <tt class="docutils literal"><span class="pre"><operators></span></tt> element, we can define any number of operators. Operators are used to define functions in the transformed space of each dimension, which in this case is Fourier space. The derivative of a function is equivalent to multiplying by <span class="math">\(i*k\)</span> in Fourier space, so the <span class="math">\(\frac{i}{2}\frac{\partial^2 \phi}{\partial \tau^2}\)</span> term in our equation of motion is equivalent to multiplying by <span class="math">\(-\frac{i}{2}k_\tau^2\)</span> in Fourier space. In this example we define “Ltt” as an operator of exactly that form, and in the equation of motion it is applied to the field “phi”.</p>
<p>Operators can be explicit (<tt class="docutils literal"><span class="pre">kind="ex"</span></tt>) or in the interaction picture (<tt class="docutils literal"><span class="pre">kind="ip"</span></tt>). The interaction picture can be more efficient, but it restricts the possible syntax of the equation of motion. Safe utilisation of interaction picture operators will be described later, but for now let us emphasise that <strong>explicit operators should be used</strong> unless the user is clear what they are doing. That said, <strong>XMDS2</strong> will generate an error if the user tries to use interaction picture operators incorrectly. The <tt class="docutils literal"><span class="pre">constant="yes"</span></tt> option in the operator block means that the operator is not a function of the propagation dimension “xi”, and therefore only needs to be calculated once at the start of the simulation.</p>
<p>The output of a partial differential equation offers more possibilities than an ordinary differential equation, and we examine some in this example.</p>
<p>For vectors with transverse dimensions, we can sample functions of the vectors on the full lattice or a subset of the points. In the <tt class="docutils literal"><span class="pre"><sampling_group></span></tt> element, we must add a string called “basis” that determines the space in which each transverse dimension is to be sampled, optionally followed by the number of points to be sampled in parentheses. If the number of points is not specified, it will default to a complete sampling of all points in that dimension. If a non-zero number of points is specified, it must be a factor of the lattice size for that dimension.</p>
<div class="highlight-xpdeint"><div class="highlight"><pre><span class="nt"><sampling_group</span> <span class="na">basis=</span><span class="s">"tau"</span> <span class="na">initial_sample=</span><span class="s">"yes"</span><span class="nt">></span>
<span class="nt"><moments></span>density<span class="nt"></moments></span>
<span class="nt"><dependencies></span>wavefunction<span class="nt"></dependencies></span>
<span class="cp"><![CDATA[</span>
<span class="n">density</span> <span class="o">=</span> <span class="nf">mod2</span><span class="p">(</span><span class="n">phi</span><span class="p">);</span>
<span class="cp">]]></span>
<span class="nt"></sampling_group></span>
</pre></div>
</div>
<p>The first output group samples the mod square of the vector “phi” over the full lattice of 128 points.</p>
<p>If the lattice parameter is set to zero points, then the corresponding dimension is integrated.</p>
<div class="highlight-xpdeint"><div class="highlight"><pre><span class="nt"><sampling_group</span> <span class="na">basis=</span><span class="s">"tau(0)"</span> <span class="na">initial_sample=</span><span class="s">"yes"</span><span class="nt">></span>
<span class="nt"><moments></span>normalisation<span class="nt"></moments></span>
<span class="nt"><dependencies></span>wavefunction<span class="nt"></dependencies></span>
<span class="cp"><![CDATA[</span>
<span class="n">normalisation</span> <span class="o">=</span> <span class="nf">mod2</span><span class="p">(</span><span class="n">phi</span><span class="p">);</span>
<span class="cp">]]></span>
<span class="nt"></sampling_group></span>
</pre></div>
</div>
<p>This second output group samples the normalisation of the wavefunction <span class="math">\(\int d\tau |\phi(\tau)|^2\)</span> over the domain of <span class="math">\(\tau\)</span>. This output requires only a single real number per sample, so in the integrate element we have chosen to sample it many more times than the vectors themselves.</p>
<p>Finally, functions of the vectors can be sampled with their dimensions in Fourier space.</p>
<div class="highlight-xpdeint"><div class="highlight"><pre><span class="nt"><sampling_group</span> <span class="na">basis=</span><span class="s">"ktau(32)"</span> <span class="na">initial_sample=</span><span class="s">"yes"</span><span class="nt">></span>
<span class="nt"><moments></span>densityK<span class="nt"></moments></span>
<span class="nt"><dependencies></span>wavefunction<span class="nt"></dependencies></span>
<span class="cp"><![CDATA[</span>
<span class="n">densityK</span> <span class="o">=</span> <span class="nf">mod2</span><span class="p">(</span><span class="n">phi</span><span class="p">);</span>
<span class="cp">]]></span>
<span class="nt"></sampling_group></span>
</pre></div>
</div>
<p>The final output group above samples the mod square of the Fourier-space wavefunction phi on a sample of 32 points.</p>
</div>
<div class="section" id="kubo-oscillator">
<span id="kubo"></span><span id="index-2"></span><h2>Kubo Oscillator<a class="headerlink" href="#kubo-oscillator" title="Permalink to this headline">¶</a></h2>
<p>This example demonstrates the integration of a stochastic differential equation. We examine the Kubo oscillator, which is a complex variable whose phase is evolving according to a Wiener noise. In a suitable rotating frame, the equation of motion for the variable is</p>
<div class="math">
\[\frac{dz}{dt} = i z \;\eta\]</div>
<p>where <span class="math">\(\eta(t)\)</span> is the Wiener differential, and we interpret this as a Stratonovich equation. In other common notation, this is sometimes written:</p>
<div class="math">
\[dz = i z \;\circ dW\]</div>
<p>Most algorithms employed by XMDS require the equations to be input in the Stratonovich form. Ito differential equations can always be transformed into Stratonovich equations, and in this case the difference is equivalent to the choice of rotating frame. This equation is solved by the following XMDS2 script:</p>
<div class="highlight-xpdeint"><div class="highlight"><pre><span class="nt"><simulation</span> <span class="na">xmds-version=</span><span class="s">"2"</span><span class="nt">></span>
<span class="nt"><name></span>kubo<span class="nt"></name></span>
<span class="nt"><author></span>Graham Dennis and Joe Hope<span class="nt"></author></span>
<span class="nt"><description></span>
Example Kubo oscillator simulation
<span class="nt"></description></span>
<span class="nt"><geometry></span>
<span class="nt"><propagation_dimension></span> t <span class="nt"></propagation_dimension></span>
<span class="nt"></geometry></span>
<span class="nt"><driver</span> <span class="na">name=</span><span class="s">"multi-path"</span> <span class="na">paths=</span><span class="s">"10000"</span> <span class="nt">/></span>
<span class="nt"><features></span>
<span class="nt"><error_check</span> <span class="nt">/></span>
<span class="nt"><benchmark</span> <span class="nt">/></span>
<span class="nt"></features></span>
<span class="nt"><noise_vector</span> <span class="na">name=</span><span class="s">"drivingNoise"</span> <span class="na">dimensions=</span><span class="s">""</span> <span class="na">kind=</span><span class="s">"wiener"</span> <span class="na">type=</span><span class="s">"real"</span> <span class="na">method=</span><span class="s">"dsfmt"</span> <span class="na">seed=</span><span class="s">"314 159 276"</span><span class="nt">></span>
<span class="nt"><components></span>eta<span class="nt"></components></span>
<span class="nt"></noise_vector></span>
<span class="nt"><vector</span> <span class="na">name=</span><span class="s">"main"</span> <span class="na">type=</span><span class="s">"complex"</span><span class="nt">></span>
<span class="nt"><components></span> z <span class="nt"></components></span>
<span class="nt"><initialisation></span>
<span class="cp"><![CDATA[</span>
<span class="n">z</span> <span class="o">=</span> <span class="mf">1.0</span><span class="p">;</span>
<span class="cp">]]></span>
<span class="nt"></initialisation></span>
<span class="nt"></vector></span>
<span class="nt"><sequence></span>
<span class="nt"><integrate</span> <span class="na">algorithm=</span><span class="s">"SI"</span> <span class="na">interval=</span><span class="s">"10"</span> <span class="na">steps=</span><span class="s">"1000"</span><span class="nt">></span>
<span class="nt"><samples></span>100<span class="nt"></samples></span>
<span class="nt"><operators></span>
<span class="nt"><integration_vectors></span>main<span class="nt"></integration_vectors></span>
<span class="nt"><dependencies></span>drivingNoise<span class="nt"></dependencies></span>
<span class="cp"><![CDATA[</span>
<span class="n">dz_dt</span> <span class="o">=</span> <span class="kc">i</span><span class="o">*</span><span class="n">z</span><span class="o">*</span><span class="n">eta</span><span class="p">;</span>
<span class="cp">]]></span>
<span class="nt"></operators></span>
<span class="nt"></integrate></span>
<span class="nt"></sequence></span>
<span class="nt"><output></span>
<span class="nt"><sampling_group</span> <span class="na">initial_sample=</span><span class="s">"yes"</span><span class="nt">></span>
<span class="nt"><moments></span>zR zI<span class="nt"></moments></span>
<span class="nt"><dependencies></span>main<span class="nt"></dependencies></span>
<span class="cp"><![CDATA[</span>
<span class="n">zR</span> <span class="o">=</span> <span class="n">z</span><span class="p">.</span><span class="n">Re</span><span class="p">();</span>
<span class="n">zI</span> <span class="o">=</span> <span class="n">z</span><span class="p">.</span><span class="n">Im</span><span class="p">();</span>
<span class="cp">]]></span>
<span class="nt"></sampling_group></span>
<span class="nt"></output></span>
<span class="nt"></simulation></span>
</pre></div>
</div>
<p id="index-3">The first new item in this script is the <tt class="docutils literal"><span class="pre"><driver></span></tt> element. This element enables us to change top level management of the simulation. Without this element, XMDS2 will integrate the stochastic equation as described. With this element and the option <tt class="docutils literal"><span class="pre">name="multi-path"</span></tt>, it will integrate it multiple times, using different random numbers each time. The output will then contain the mean values and standard errors of your output variables. The number of integrations included in the averages is set with the <tt class="docutils literal"><span class="pre">paths</span></tt> variable.</p>
<p>In the <tt class="docutils literal"><span class="pre"><features></span></tt> element we have included the <tt class="docutils literal"><span class="pre"><error_check></span></tt> element. This performs the integration first with the specified number of steps (or with the specified tolerance), and then with twice the number of steps (or equivalently reduced tolerance). The output then includes the difference between the output variables on the coarse and the fine grids as the ‘error’ in the output variables. This error is particularly useful for stochastic integrations, where algorithms with adaptive step-sizes are less safe, so the number of integration steps must be user-specified.</p>
<p id="index-4">We define the stochastic elements in a simulation with the <tt class="docutils literal"><span class="pre"><noise_vector></span></tt> element.</p>
<div class="highlight-xpdeint"><div class="highlight"><pre><span class="nt"><noise_vector</span> <span class="na">name=</span><span class="s">"drivingNoise"</span> <span class="na">dimensions=</span><span class="s">""</span> <span class="na">kind=</span><span class="s">"wiener"</span> <span class="na">type=</span><span class="s">"real"</span> <span class="na">method=</span><span class="s">"dsfmt"</span> <span class="na">seed=</span><span class="s">"314 159 276"</span><span class="nt">></span>
<span class="nt"><components></span>eta<span class="nt"></components></span>
<span class="nt"></noise_vector></span>
</pre></div>
</div>
<p>This defines a vector that is used like any other, but it will be randomly generated with particular statistics and characteristics rather than initialised. The name, dimensions and type tags are defined just as for normal vectors. The names of the components are also defined in the same way. The noise is defined as a Wiener noise here (<tt class="docutils literal"><span class="pre">kind</span> <span class="pre">=</span> <span class="pre">"wiener"</span></tt>), which is a zero-mean Gaussian random noise with an average variance equal to the discretisation volume (here it is just the step size in the propagation dimension, as it is not defined over transverse dimensions). Other noise types are possible, including uniform and Poissonian noises, but we will not describe them in detail here.</p>
<p>We may also define a noise method to choose a non-default pseudo random number generator, and a seed for the random number generator. Using a seed can be very useful when debugging the behaviour of a simulation, and many compilers have pseudo-random number generators that are superior to the default option (posix).</p>
<p>The integrate block is using the semi-implicit algorithm (<tt class="docutils literal"><span class="pre">algorithm="SI"</span></tt>), which is a good default choice for stochastic problems, even though it is only second order convergent for deterministic equations. More will be said about algorithm choice later, but for now we should note that adaptive algorithms based on Runge-Kutta methods are not guaranteed to converge safely for stochastic equations. This can be particularly deceptive as they often succeed, particularly for almost any problem for which there is a known analytic solution.</p>
<p>We include elements from the noise vector in the equation of motion just as we do for any other vector. The default SI and Runge-Kutta algorithms converge to the <em>Stratonovich</em> integral. Ito stochastic equations can be converted to Stratonovich form and vice versa.</p>
<p>Executing the generated program ‘kubo’ gives slightly different output due to the “multi-path” driver.</p>
<div class="highlight-none"><div class="highlight"><pre>$ ./kubo
Beginning full step integration ...
Starting path 1
Starting path 2
... many lines omitted ...
Starting path 9999
Starting path 10000
Beginning half step integration ...
Starting path 1
Starting path 2
... many lines omitted ...
Starting path 9999
Starting path 10000
Generating output for kubo
Maximum step error in moment group 1 was 4.942549e-04
Time elapsed for simulation is: 2.71 seconds
</pre></div>
</div>
<p>The maximum step error in each moment group is given in absolute terms. This is the largest difference between the full step integration and the half step integration. While a single path might be very stochastic:</p>
<div class="figure align-center">
<img alt="_images/kuboSingle.png" src="_images/kuboSingle.png" />
<p class="caption">The mean value of the real and imaginary components of the z variable for a single path of the simulation.</p>
</div>
<p>The average over multiple paths can be increasingly smooth.</p>
<div class="figure align-center">
<img alt="_images/kubo10000.png" src="_images/kubo10000.png" />
<p class="caption">The mean and standard error of the z variable averaged over 10000 paths, as given by this simulation. It agrees within the standard error with the expected result of <span class="math">\(\exp(-t/2)\)</span>.</p>
</div>
</div>
<div class="section" id="fibre-noise">
<span id="fibre"></span><span id="index-5"></span><h2>Fibre Noise<a class="headerlink" href="#fibre-noise" title="Permalink to this headline">¶</a></h2>
<p>This simulation is a stochastic partial differential equation, in which a one-dimensional damped field is subject to a complex noise. This script can be found in <tt class="docutils literal"><span class="pre">examples/fibre.xmds</span></tt>.</p>
<div class="math">
\[\frac{\partial \psi}{\partial t} = -i \frac{\partial^2 \psi}{\partial x^2} -\gamma \psi+\beta \frac{1}{\sqrt{2}}\left(\eta_1(x)+i\eta_2(x)\right)\]</div>
<p>where the noise terms <span class="math">\(\eta_j(x,t)\)</span> are Wiener differentials and the equation is interpreted as a Stratonovich differential equation. On a finite grid, these increments have variance <span class="math">\(\frac{1}{\Delta x \Delta t}\)</span>.</p>
<div class="highlight-xpdeint"><div class="highlight"><pre><span class="nt"><simulation</span> <span class="na">xmds-version=</span><span class="s">"2"</span><span class="nt">></span>
<span class="nt"><name></span>fibre<span class="nt"></name></span>
<span class="nt"><author></span>Joe Hope and Graham Dennis<span class="nt"></author></span>
<span class="nt"><description></span>
Example fibre noise simulation
<span class="nt"></description></span>
<span class="nt"><geometry></span>
<span class="nt"><propagation_dimension></span> t <span class="nt"></propagation_dimension></span>
<span class="nt"><transverse_dimensions></span>
<span class="nt"><dimension</span> <span class="na">name=</span><span class="s">"x"</span> <span class="na">lattice=</span><span class="s">"64"</span> <span class="na">domain=</span><span class="s">"(-5, 5)"</span> <span class="nt">/></span>
<span class="nt"></transverse_dimensions></span>
<span class="nt"></geometry></span>
<span class="nt"><driver</span> <span class="na">name=</span><span class="s">"mpi-multi-path"</span> <span class="na">paths=</span><span class="s">"8"</span> <span class="nt">/></span>
<span class="nt"><features></span>
<span class="nt"><auto_vectorise</span> <span class="nt">/></span>
<span class="nt"><benchmark</span> <span class="nt">/></span>
<span class="nt"><error_check</span> <span class="nt">/></span>
<span class="nt"><globals></span>
<span class="cp"><![CDATA[</span>
<span class="k">const</span> <span class="nf">real</span> <span class="n">ggamma</span> <span class="o">=</span> <span class="mf">1.0</span><span class="p">;</span>
<span class="k">const</span> <span class="nf">real</span> <span class="n">beta</span> <span class="o">=</span> <span class="n">sqrt</span><span class="p">(</span><span class="n">M_PI</span><span class="o">*</span><span class="n">ggamma</span><span class="o">/</span><span class="mf">10.0</span><span class="p">);</span>
<span class="cp">]]></span>
<span class="nt"></globals></span>
<span class="nt"></features></span>
<span class="nt"><noise_vector</span> <span class="na">name=</span><span class="s">"drivingNoise"</span> <span class="na">dimensions=</span><span class="s">"x"</span> <span class="na">kind=</span><span class="s">"wiener"</span> <span class="na">type=</span><span class="s">"complex"</span> <span class="na">method=</span><span class="s">"dsfmt"</span> <span class="na">seed=</span><span class="s">"314 159 276"</span><span class="nt">></span>
<span class="nt"><components></span>Eta<span class="nt"></components></span>
<span class="nt"></noise_vector></span>
<span class="nt"><vector</span> <span class="na">name=</span><span class="s">"main"</span> <span class="na">initial_basis=</span><span class="s">"x"</span> <span class="na">type=</span><span class="s">"complex"</span><span class="nt">></span>
<span class="nt"><components></span>phi<span class="nt"></components></span>
<span class="nt"><initialisation></span>
<span class="cp"><![CDATA[</span>
<span class="n">phi</span> <span class="o">=</span> <span class="mf">0.0</span><span class="p">;</span>
<span class="cp">]]></span>
<span class="nt"></initialisation></span>
<span class="nt"></vector></span>
<span class="nt"><sequence></span>
<span class="nt"><integrate</span> <span class="na">algorithm=</span><span class="s">"SI"</span> <span class="na">iterations=</span><span class="s">"3"</span> <span class="na">interval=</span><span class="s">"2.5"</span> <span class="na">steps=</span><span class="s">"200000"</span><span class="nt">></span>
<span class="nt"><samples></span>50<span class="nt"></samples></span>
<span class="nt"><operators></span>
<span class="nt"><operator</span> <span class="na">kind=</span><span class="s">"ex"</span><span class="nt">></span>
<span class="nt"><operator_names></span>L<span class="nt"></operator_names></span>
<span class="cp"><![CDATA[</span>
<span class="n">L</span> <span class="o">=</span> <span class="o">-</span><span class="kc">i</span><span class="o">*</span><span class="n">kx</span><span class="o">*</span><span class="n">kx</span><span class="p">;</span>
<span class="cp">]]></span>
<span class="nt"></operator></span>
<span class="nt"><dependencies></span>drivingNoise<span class="nt"></dependencies></span>
<span class="nt"><integration_vectors></span>main<span class="nt"></integration_vectors></span>
<span class="cp"><![CDATA[</span>
<span class="n">dphi_dt</span> <span class="o">=</span> <span class="n">L</span><span class="p">[</span><span class="n">phi</span><span class="p">]</span> <span class="o">-</span> <span class="n">ggamma</span><span class="o">*</span><span class="n">phi</span> <span class="o">+</span> <span class="n">beta</span><span class="o">*</span><span class="n">Eta</span><span class="p">;</span>
<span class="cp">]]></span>
<span class="nt"></operators></span>
<span class="nt"></integrate></span>
<span class="nt"></sequence></span>
<span class="nt"><output></span>
<span class="nt"><sampling_group</span> <span class="na">basis=</span><span class="s">"kx"</span> <span class="na">initial_sample=</span><span class="s">"yes"</span><span class="nt">></span>
<span class="nt"><moments></span>pow_dens<span class="nt"></moments></span>
<span class="nt"><dependencies></span>main<span class="nt"></dependencies></span>
<span class="cp"><![CDATA[</span>
<span class="n">pow_dens</span> <span class="o">=</span> <span class="nf">mod2</span><span class="p">(</span><span class="n">phi</span><span class="p">);</span>
<span class="cp">]]></span>
<span class="nt"></sampling_group></span>
<span class="nt"></output></span>
<span class="nt"></simulation></span>
</pre></div>
</div>
<p>Note that the noise vector used in this example is complex-valued, and has the argument <tt class="docutils literal"><span class="pre">dimensions="x"</span></tt> to define it as a field of delta-correlated noises along the x-dimension.</p>
<span class="target" id="index-6"></span><p id="index-7">This simulation demonstrates the ease with which XMDS2 can be used in a parallel processing environment. Instead of using the stochastic driver “multi-path”, we simply replace it with “mpi-multi-path”. This instructs XMDS2 to write a parallel version of the program based on the widespread <a class="reference external" href="http://www.open-mpi.org/">MPI standard</a>. This protocol allows multiple processors or clusters of computers to work simultaneously on the same problem. Free open source libraries implementing this standard can be installed on a linux machine, and come standard on Mac OS X. They are also common on many supercomputer architectures. Parallel processing can also be used with deterministic problems to great effect, as discussed in the later example <a class="reference internal" href="#wignerarguments"><em>Wigner Function</em></a>.</p>
<p>Executing this program is slightly different with the MPI option. The details can change between MPI implementations, but as an example:</p>
<div class="highlight-none"><div class="highlight"><pre>$xmds2 fibre.xmds
xmds2 version 2.1 "Happy Mollusc" (r2543)
Copyright 2000-2012 Graham Dennis, Joseph Hope, Mattias Johnsson
and the xmds team
Generating source code...
... done
Compiling simulation...
... done. Type './fibre' to run.
</pre></div>
</div>
<p>Note that different compile options (and potentially a different compiler) are used by XMDS2, but this is transparent to the user. MPI simulations will have to be run using syntax that will depend on the MPI implementation. Here we show the version based on the popular open source <a class="reference external" href="http://www.open-mpi.org/">Open-MPI</a> implementation.</p>
<div class="highlight-none"><div class="highlight"><pre>$ mpirun -np 4 ./fibre
Found enlightenment... (Importing wisdom)
Planning for x <---> kx transform... done.
Beginning full step integration ...
Rank[0]: Starting path 1
Rank[1]: Starting path 2
Rank[2]: Starting path 3
Rank[3]: Starting path 4
Rank[3]: Starting path 8
Rank[0]: Starting path 5
Rank[1]: Starting path 6
Rank[2]: Starting path 7
Rank[3]: Starting path 4
Beginning half step integration ...
Rank[0]: Starting path 1
Rank[2]: Starting path 3
Rank[1]: Starting path 2
Rank[3]: Starting path 8
Rank[0]: Starting path 5
Rank[2]: Starting path 7
Rank[1]: Starting path 6
Generating output for fibre
Maximum step error in moment group 1 was 4.893437e-04
Time elapsed for simulation is: 20.99 seconds
</pre></div>
</div>
<p>In this example we used four processors. The different processors are labelled by their “Rank”, starting at zero. Because the processors are working independently, the output from the different processors can come in a randomised order. In the end, however, the .xsil and data files are constructed identically to the single processor outputs.</p>
<p>The analytic solution to the stochastic averages of this equation is given by</p>
<div class="math">
\[\langle |\psi(k,t)|^2 \rangle = \exp(-2\gamma t)|\psi(k,0)|^2 +\frac{\beta^2 L_x}{4\pi \gamma} \left(1-\exp(-2\gamma t)\right)\]</div>
<p>where <span class="math">\(L_x\)</span> is the length of the x domain. We see that a single integration of these equations is quite chaotic:</p>
<div class="figure align-center">
<img alt="_images/fibreSingle.png" src="_images/fibreSingle.png" />
<p class="caption">The momentum space density of the field as a function of time for a single path realisation.</p>
</div>
<p>while an average of 1024 paths (change <tt class="docutils literal"><span class="pre">paths="8"</span></tt> to <tt class="docutils literal"><span class="pre">paths="1024"</span></tt> in the <tt class="docutils literal"><span class="pre"><driver></span></tt> element) converges nicely to the analytic solution:</p>
<div class="figure align-center">
<img alt="_images/fibre1024.png" src="_images/fibre1024.png" />
<p class="caption">The momentum space density of the field as a function of time for an average of 1024 paths.</p>
</div>
</div>
<div class="section" id="integer-dimensions">
<span id="integerdimensionexample"></span><span id="index-8"></span><h2>Integer Dimensions<a class="headerlink" href="#integer-dimensions" title="Permalink to this headline">¶</a></h2>
<p>This example shows how to handle systems with integer-valued transverse dimensions. We will integrate the following set of equations</p>
<div class="math">
\[\frac{dx_j}{dt} = x_j \left(x_{j-1}-x_{j+1}\right)\]</div>
<p>where <span class="math">\(x_j\)</span> are complex-valued variables defined on a ring, such that <span class="math">\(j\in \{0,j_{max}\}\)</span> and the <span class="math">\(x_{j_{max}+1}\)</span> variable is identified with the variable <span class="math">\(x_{0}\)</span>, and the variable <span class="math">\(x_{-1}\)</span> is identified with the variable <span class="math">\(x_{j_{max}}\)</span>.</p>
<div class="highlight-xpdeint"><div class="highlight"><pre><span class="nt"><simulation</span> <span class="na">xmds-version=</span><span class="s">"2"</span><span class="nt">></span>
<span class="nt"><name></span>integer_dimensions<span class="nt"></name></span>
<span class="nt"><author></span>Graham Dennis<span class="nt"></author></span>
<span class="nt"><description></span>
XMDS2 script to test integer dimensions.
<span class="nt"></description></span>
<span class="nt"><features></span>
<span class="nt"><benchmark</span> <span class="nt">/></span>
<span class="nt"><error_check</span> <span class="nt">/></span>
<span class="nt"><bing</span> <span class="nt">/></span>
<span class="nt"><diagnostics</span> <span class="nt">/></span> <span class="c"><!-- This will make sure that all nonlocal accesses of dimensions are safe --></span>
<span class="nt"></features></span>
<span class="nt"><geometry></span>
<span class="nt"><propagation_dimension></span> t <span class="nt"></propagation_dimension></span>
<span class="nt"><transverse_dimensions></span>
<span class="nt"><dimension</span> <span class="na">name=</span><span class="s">"j"</span> <span class="na">type=</span><span class="s">"integer"</span> <span class="na">lattice=</span><span class="s">"5"</span> <span class="na">domain=</span><span class="s">"(0,4)"</span> <span class="nt">/></span>
<span class="nt"></transverse_dimensions></span>
<span class="nt"></geometry></span>
<span class="nt"><vector</span> <span class="na">name=</span><span class="s">"main"</span> <span class="na">type=</span><span class="s">"complex"</span><span class="nt">></span>
<span class="nt"><components></span> x <span class="nt"></components></span>
<span class="nt"><initialisation></span>
<span class="cp"><![CDATA[</span>
<span class="n">x</span> <span class="o">=</span> <span class="mf">1.0e-3</span><span class="p">;</span>
<span class="n">x</span><span class="p">(</span><span class="n">j</span> <span class="o">=></span> <span class="mi">0</span><span class="p">)</span> <span class="o">=</span> <span class="mf">1.0</span><span class="p">;</span>
<span class="cp">]]></span>
<span class="nt"></initialisation></span>
<span class="nt"></vector></span>
<span class="nt"><sequence></span>
<span class="nt"><integrate</span> <span class="na">algorithm=</span><span class="s">"ARK45"</span> <span class="na">interval=</span><span class="s">"60"</span> <span class="na">steps=</span><span class="s">"25000"</span> <span class="na">tolerance=</span><span class="s">"1.0e-9"</span><span class="nt">></span>
<span class="nt"><samples></span>1000<span class="nt"></samples></span>
<span class="nt"><operators></span>
<span class="nt"><integration_vectors></span>main<span class="nt"></integration_vectors></span>
<span class="cp"><![CDATA[</span>
<span class="kt">long</span> <span class="n">j_minus_one</span> <span class="o">=</span> <span class="p">(</span><span class="n">j</span><span class="o">-</span><span class="mi">1</span><span class="p">)</span> <span class="o">%</span> <span class="n">_lattice_j</span><span class="p">;</span>
<span class="k">if</span> <span class="p">(</span><span class="n">j_minus_one</span> <span class="o"><</span> <span class="mi">0</span><span class="p">)</span>
<span class="n">j_minus_one</span> <span class="o">+=</span> <span class="n">_lattice_j</span><span class="p">;</span>
<span class="kt">long</span> <span class="n">j_plus_one</span> <span class="o">=</span> <span class="p">(</span><span class="n">j</span><span class="o">+</span><span class="mi">1</span><span class="p">)</span> <span class="o">%</span> <span class="n">_lattice_j</span><span class="p">;</span>
<span class="n">dx_dt</span><span class="p">(</span><span class="n">j</span> <span class="o">=></span> <span class="n">j</span><span class="p">)</span> <span class="o">=</span> <span class="n">x</span><span class="p">(</span><span class="n">j</span> <span class="o">=></span> <span class="n">j</span><span class="p">)</span><span class="o">*</span><span class="p">(</span><span class="n">x</span><span class="p">(</span><span class="n">j</span> <span class="o">=></span> <span class="n">j_minus_one</span><span class="p">)</span> <span class="o">-</span> <span class="n">x</span><span class="p">(</span><span class="n">j</span> <span class="o">=></span> <span class="n">j_plus_one</span><span class="p">));</span>
<span class="cp">]]></span>
<span class="nt"></operators></span>
<span class="nt"></integrate></span>
<span class="nt"></sequence></span>
<span class="nt"><output></span>
<span class="nt"><sampling_group</span> <span class="na">basis=</span><span class="s">"j"</span> <span class="na">initial_sample=</span><span class="s">"yes"</span><span class="nt">></span>
<span class="nt"><moments></span>xR<span class="nt"></moments></span>
<span class="nt"><dependencies></span>main<span class="nt"></dependencies></span>
<span class="cp"><![CDATA[</span>
<span class="n">xR</span> <span class="o">=</span> <span class="n">x</span><span class="p">.</span><span class="n">Re</span><span class="p">();</span>
<span class="cp">]]></span>
<span class="nt"></sampling_group></span>
<span class="nt"></output></span>
<span class="nt"></simulation></span>
</pre></div>
</div>
<p>The first extra feature we have used in this script is the <tt class="docutils literal"><span class="pre"><diagnostics></span></tt> element. It performs run-time checking that our generated code does not accidentally attempt to access a part of our vector that does not exist. Removing this tag will increase the speed of the simulation, but its presence helps catch coding errors.</p>
<p>The simulation defines a vector with a single transverse dimension labelled “j”, of type “integer” (“int” and “long” can also be used as synonyms for “integer”). In the absence of an explicit type, the dimension is assumed to be real-valued. The dimension has a “domain” argument as normal, defining the minimum and maximum values of the dimension’s range. The lattice element, if specified, is used as a check on the size of the domain, and will create an error if the two do not match.</p>
<p>Integer-valued dimensions can be called non-locally. Real-valued dimensions are typically coupled non-locally only through local operations in the transformed space of the dimension, but can be called non-locally in certain other situations as described in <a class="reference internal" href="reference_elements.html#referencingnonlocal"><em>the reference</em></a>. The syntax for calling integer dimensions non-locally can be seen in the initialisation CDATA block:</p>
<div class="highlight-xpdeint"><div class="highlight"><pre>x = 1.0e-3;
x(j => 0) = 1.0;
</pre></div>
</div>
<p>where the syntax <tt class="docutils literal"><span class="pre">x(j</span> <span class="pre">=></span> <span class="pre">0)</span></tt> is used to reference the variable <span class="math">\(x_0\)</span> directly. We see a more elaborate example in the integrate CDATA block:</p>
<div class="highlight-xpdeint"><div class="highlight"><pre>dx_dt(j => j) = x(j => j)*(x(j => j_minus_one) - x(j => j_plus_one));
</pre></div>
</div>
<p>where the vector “x” is called using locally defined variables. This syntax is chosen so that multiple dimensions can be addressed non-locally with minimal possibility for confusion.</p>
</div>
<div class="section" id="wigner-function">
<span id="wignerarguments"></span><h2>Wigner Function<a class="headerlink" href="#wigner-function" title="Permalink to this headline">¶</a></h2>
<p>This example integrates the two-dimensional partial differential equation</p>
<div class="math">
\[\begin{split}\begin{split}
\frac{\partial W}{\partial t} &= \Bigg[ \left(\omega + \frac{U_{int}}{\hbar}\left(x^2+y^2-1\right)\right) \left(x \frac{\partial}{\partial y}
- y \frac{\partial}{\partial x}\right)\\
&\phantom{=\Bigg[} - \frac{U_{int}}{16 \hbar}\left(x\left(\frac{\partial^3}{\partial x^2 \partial y}
+\frac{\partial^3}{\partial y^3}\right)-y\left(\frac{\partial^3}{\partial y^2 \partial x}+\frac{\partial^3}{\partial x^3}\right)\right)\Bigg]W(x,y,t)
\end{split}\end{split}\]</div>
<p>with the added restriction that the derivative is forced to zero outside a certain radius. This extra condition helps maintain the long-term stability of the integration. The script can be found in <tt class="docutils literal"><span class="pre">examples/wigner_arguments_mpi.xmds</span></tt> under your XMDS2 installation directory.</p>
<div class="highlight-xpdeint"><div class="highlight"><pre><span class="nt"><simulation</span> <span class="na">xmds-version=</span><span class="s">"2"</span><span class="nt">></span>
<span class="nt"><name></span>wigner<span class="nt"></name></span>
<span class="nt"><author></span>Graham Dennis and Joe Hope<span class="nt"></author></span>
<span class="nt"><description></span>
Simulation of the Wigner function for an anharmonic oscillator with the initial state
being a coherent state.
<span class="nt"></description></span>
<span class="nt"><features></span>
<span class="nt"><benchmark</span> <span class="nt">/></span>
<span class="nt"><globals></span>
<span class="cp"><![CDATA[</span>
<span class="nf">real</span> <span class="n">Uint_hbar_on16</span><span class="p">;</span>
<span class="cp">]]></span>
<span class="nt"></globals></span>
<span class="nt"><arguments></span>
<span class="nt"><argument</span> <span class="na">name=</span><span class="s">"omega"</span> <span class="na">type=</span><span class="s">"real"</span> <span class="na">default_value=</span><span class="s">"0.0"</span> <span class="nt">/></span>
<span class="nt"><argument</span> <span class="na">name=</span><span class="s">"alpha_0"</span> <span class="na">type=</span><span class="s">"real"</span> <span class="na">default_value=</span><span class="s">"3.0"</span> <span class="nt">/></span>
<span class="nt"><argument</span> <span class="na">name=</span><span class="s">"absorb"</span> <span class="na">type=</span><span class="s">"real"</span> <span class="na">default_value=</span><span class="s">"8.0"</span> <span class="nt">/></span>
<span class="nt"><argument</span> <span class="na">name=</span><span class="s">"width"</span> <span class="na">type=</span><span class="s">"real"</span> <span class="na">default_value=</span><span class="s">"0.3"</span> <span class="nt">/></span>
<span class="nt"><argument</span> <span class="na">name=</span><span class="s">"Uint_hbar"</span> <span class="na">type=</span><span class="s">"real"</span> <span class="na">default_value=</span><span class="s">"1.0"</span> <span class="nt">/></span>
<span class="cp"><![CDATA[</span>
<span class="cm">/* derived constants */</span>
<span class="n">Uint_hbar_on16</span> <span class="o">=</span> <span class="n">Uint_hbar</span><span class="o">/</span><span class="mf">16.0</span><span class="p">;</span>
<span class="cp">]]></span>
<span class="nt"></arguments></span>
<span class="nt"><bing</span> <span class="nt">/></span>
<span class="nt"><fftw</span> <span class="na">plan=</span><span class="s">"patient"</span> <span class="nt">/></span>
<span class="nt"><openmp</span> <span class="nt">/></span>
<span class="nt"></features></span>
<span class="nt"><driver</span> <span class="na">name=</span><span class="s">"distributed-mpi"</span> <span class="nt">/></span>
<span class="nt"><geometry></span>
<span class="nt"><propagation_dimension></span> t <span class="nt"></propagation_dimension></span>
<span class="nt"><transverse_dimensions></span>
<span class="nt"><dimension</span> <span class="na">name=</span><span class="s">"x"</span> <span class="na">lattice=</span><span class="s">"128"</span> <span class="na">domain=</span><span class="s">"(-6, 6)"</span> <span class="nt">/></span>
<span class="nt"><dimension</span> <span class="na">name=</span><span class="s">"y"</span> <span class="na">lattice=</span><span class="s">"128"</span> <span class="na">domain=</span><span class="s">"(-6, 6)"</span> <span class="nt">/></span>
<span class="nt"></transverse_dimensions></span>
<span class="nt"></geometry></span>
<span class="nt"><vector</span> <span class="na">name=</span><span class="s">"main"</span> <span class="na">initial_basis=</span><span class="s">"x y"</span> <span class="na">type=</span><span class="s">"complex"</span><span class="nt">></span>
<span class="nt"><components></span> W <span class="nt"></components></span>
<span class="nt"><initialisation></span>
<span class="cp"><![CDATA[</span>
<span class="n">W</span> <span class="o">=</span> <span class="mf">2.0</span><span class="o">/</span><span class="n">M_PI</span> <span class="o">*</span> <span class="n">exp</span><span class="p">(</span><span class="o">-</span><span class="mf">2.0</span><span class="o">*</span><span class="p">(</span><span class="n">y</span><span class="o">*</span><span class="n">y</span> <span class="o">+</span> <span class="p">(</span><span class="n">x</span><span class="o">-</span><span class="n">alpha_0</span><span class="p">)</span><span class="o">*</span><span class="p">(</span><span class="n">x</span><span class="o">-</span><span class="n">alpha_0</span><span class="p">)));</span>
<span class="cp">]]></span>
<span class="nt"></initialisation></span>
<span class="nt"></vector></span>
<span class="nt"><vector</span> <span class="na">name=</span><span class="s">"dampConstants"</span> <span class="na">initial_basis=</span><span class="s">"x y"</span> <span class="na">type=</span><span class="s">"real"</span><span class="nt">></span>
<span class="nt"><components></span>damping<span class="nt"></components></span>
<span class="nt"><initialisation></span>
<span class="cp"><![CDATA[</span>
<span class="k">if</span> <span class="p">(</span><span class="n">sqrt</span><span class="p">(</span><span class="n">x</span><span class="o">*</span><span class="n">x</span> <span class="o">+</span> <span class="n">y</span><span class="o">*</span><span class="n">y</span><span class="p">)</span> <span class="o">></span> <span class="n">_max_x</span><span class="o">-</span><span class="n">width</span><span class="p">)</span>
<span class="n">damping</span> <span class="o">=</span> <span class="mf">0.0</span><span class="p">;</span>
<span class="k">else</span>
<span class="n">damping</span> <span class="o">=</span> <span class="mf">1.0</span><span class="p">;</span>
<span class="cp">]]></span>
<span class="nt"></initialisation></span>
<span class="nt"></vector></span>
<span class="nt"><sequence></span>
<span class="nt"><integrate</span> <span class="na">algorithm=</span><span class="s">"ARK89"</span> <span class="na">tolerance=</span><span class="s">"1e-7"</span> <span class="na">interval=</span><span class="s">"7.0e-4"</span> <span class="na">steps=</span><span class="s">"100000"</span><span class="nt">></span>
<span class="nt"><samples></span>50<span class="nt"></samples></span>
<span class="nt"><operators></span>
<span class="nt"><operator</span> <span class="na">kind=</span><span class="s">"ex"</span><span class="nt">></span>
<span class="nt"><operator_names></span>Lx Ly Lxxx Lxxy Lxyy Lyyy<span class="nt"></operator_names></span>
<span class="cp"><![CDATA[</span>
<span class="n">Lx</span> <span class="o">=</span> <span class="kc">i</span><span class="o">*</span><span class="n">kx</span><span class="p">;</span>
<span class="n">Ly</span> <span class="o">=</span> <span class="kc">i</span><span class="o">*</span><span class="n">ky</span><span class="p">;</span>
<span class="n">Lxxx</span> <span class="o">=</span> <span class="o">-</span><span class="kc">i</span><span class="o">*</span><span class="n">kx</span><span class="o">*</span><span class="n">kx</span><span class="o">*</span><span class="n">kx</span><span class="p">;</span>
<span class="n">Lxxy</span> <span class="o">=</span> <span class="o">-</span><span class="kc">i</span><span class="o">*</span><span class="n">kx</span><span class="o">*</span><span class="n">kx</span><span class="o">*</span><span class="n">ky</span><span class="p">;</span>
<span class="n">Lxyy</span> <span class="o">=</span> <span class="o">-</span><span class="kc">i</span><span class="o">*</span><span class="n">kx</span><span class="o">*</span><span class="n">ky</span><span class="o">*</span><span class="n">ky</span><span class="p">;</span>
<span class="n">Lyyy</span> <span class="o">=</span> <span class="o">-</span><span class="kc">i</span><span class="o">*</span><span class="n">ky</span><span class="o">*</span><span class="n">ky</span><span class="o">*</span><span class="n">ky</span><span class="p">;</span>
<span class="cp">]]></span>
<span class="nt"></operator></span>
<span class="nt"><integration_vectors></span>main<span class="nt"></integration_vectors></span>
<span class="nt"><dependencies></span>dampConstants<span class="nt"></dependencies></span>
<span class="cp"><![CDATA[</span>
<span class="nf">real</span> <span class="n">rotation</span> <span class="o">=</span> <span class="n">omega</span> <span class="o">+</span> <span class="n">Uint_hbar</span><span class="o">*</span><span class="p">(</span><span class="o">-</span><span class="mf">1.0</span> <span class="o">+</span> <span class="n">x</span><span class="o">*</span><span class="n">x</span> <span class="o">+</span> <span class="n">y</span><span class="o">*</span><span class="n">y</span><span class="p">);</span>
<span class="n">dW_dt</span> <span class="o">=</span> <span class="n">damping</span> <span class="o">*</span> <span class="p">(</span> <span class="n">rotation</span> <span class="o">*</span> <span class="p">(</span><span class="n">x</span><span class="o">*</span><span class="n">Ly</span><span class="p">[</span><span class="n">W</span><span class="p">]</span> <span class="o">-</span> <span class="n">y</span><span class="o">*</span><span class="n">Lx</span><span class="p">[</span><span class="n">W</span><span class="p">])</span>
<span class="o">-</span> <span class="n">Uint_hbar_on16</span><span class="o">*</span><span class="p">(</span> <span class="n">x</span><span class="o">*</span><span class="p">(</span><span class="n">Lxxy</span><span class="p">[</span><span class="n">W</span><span class="p">]</span> <span class="o">+</span> <span class="n">Lyyy</span><span class="p">[</span><span class="n">W</span><span class="p">])</span> <span class="o">-</span> <span class="n">y</span><span class="o">*</span><span class="p">(</span><span class="n">Lxyy</span><span class="p">[</span><span class="n">W</span><span class="p">]</span> <span class="o">+</span> <span class="n">Lxxx</span><span class="p">[</span><span class="n">W</span><span class="p">])</span> <span class="p">)</span>
<span class="p">);</span>
<span class="cp">]]></span>
<span class="nt"></operators></span>
<span class="nt"></integrate></span>
<span class="nt"></sequence></span>
<span class="nt"><output></span>
<span class="nt"><sampling_group</span> <span class="na">basis=</span><span class="s">"x y"</span> <span class="na">initial_sample=</span><span class="s">"yes"</span><span class="nt">></span>
<span class="nt"><moments></span>WR WI<span class="nt"></moments></span>
<span class="nt"><dependencies></span>main<span class="nt"></dependencies></span>
<span class="cp"><![CDATA[</span>
<span class="n">_SAMPLE_COMPLEX</span><span class="p">(</span><span class="n">W</span><span class="p">);</span>
<span class="cp">]]></span>
<span class="nt"></sampling_group></span>
<span class="nt"></output></span>
<span class="nt"></simulation></span>
</pre></div>
</div>
<span class="target" id="index-9"></span><p id="index-10">This example demonstrates two new features of XMDS2. The first is the use of parallel processing for a deterministic problem. The FFTW library only allows MPI processing of multidimensional vectors. For multidimensional simulations, the generated program can be parallelised simply by adding the <tt class="docutils literal"><span class="pre">name="distributed-mpi"</span></tt> argument to the <tt class="docutils literal"><span class="pre"><driver></span></tt> element.</p>
<div class="highlight-xpdeint"><div class="highlight"><pre>$ xmds2 wigner_argument_mpi.xmds
xmds2 version 2.1 "Happy Mollusc" (r2680)
Copyright 2000-2012 Graham Dennis, Joseph Hope, Mattias Johnsson
and the xmds team
Generating source code...
... done
Compiling simulation...
... done. Type './wigner' to run.
</pre></div>
</div>
<p>To use multiple processors, the final program is then called using the (implementation specific) MPI wrapper:</p>
<div class="highlight-xpdeint"><div class="highlight"><pre>$ mpirun -np 2 ./wigner
Planning for (distributed x, y) <span class="nt"><---></span> (distributed ky, kx) transform... done.
Planning for (distributed x, y) <span class="nt"><---></span> (distributed ky, kx) transform... done.
Sampled field (for moment group #1) at t = 0.000000e+00
Current timestep: 5.908361e-06
Sampled field (for moment group #1) at t = 1.400000e-05
Current timestep: 4.543131e-06
...
</pre></div>
</div>
<p>The possible acceleration achievable when parallelising a given simulation depends on a great many things including available memory and cache. As a general rule, it will improve as the simulation size gets larger, but the easiest way to find out is to test. The optimum speed up is obviously proportional to the number of available processing cores.</p>
<p id="index-11">The second new feature in this simulation is the <tt class="docutils literal"><span class="pre"><arguments></span></tt> element in the <tt class="docutils literal"><span class="pre"><features></span></tt> block. This is a way of specifying global variables with a given type that can then be input at run time. The variables are specified in a self explanatory way</p>
<div class="highlight-xpdeint"><div class="highlight"><pre><span class="nt"><arguments></span>
<span class="nt"><argument</span> <span class="na">name=</span><span class="s">"omega"</span> <span class="na">type=</span><span class="s">"real"</span> <span class="na">default_value=</span><span class="s">"0.0"</span> <span class="nt">/></span>
...
<span class="nt"><argument</span> <span class="na">name=</span><span class="s">"Uint_hbar"</span> <span class="na">type=</span><span class="s">"real"</span> <span class="na">default_value=</span><span class="s">"1.0"</span> <span class="nt">/></span>
<span class="nt"></arguments></span>
</pre></div>
</div>
<p>where the “default_value” is used as the valuable of the variable if no arguments are given. In the absence of the generating script, the program can document its options with the <tt class="docutils literal"><span class="pre">--help</span></tt> argument:</p>
<div class="highlight-none"><div class="highlight"><pre>$ ./wigner --help
Usage: wigner --omega <real> --alpha_0 <real> --absorb <real> --width <real> --Uint_hbar <real>
Details:
Option Type Default value
-o, --omega real 0.0
-a, --alpha_0 real 3.0
-b, --absorb real 8.0
-w, --width real 0.3
-U, --Uint_hbar real 1.0
</pre></div>
</div>
<p>We can change one or more of these variables’ values in the simulation by passing it at run time.</p>
<div class="highlight-none"><div class="highlight"><pre>$ mpirun -np 2 ./wigner --omega 0.1 --alpha_0 2.5 --Uint_hbar 0
Found enlightenment... (Importing wisdom)
Planning for (distributed x, y) <---> (distributed ky, kx) transform... done.
Planning for (distributed x, y) <---> (distributed ky, kx) transform... done.
Sampled field (for moment group #1) at t = 0.000000e+00
Current timestep: 1.916945e-04
...
</pre></div>
</div>
<p>The values that were used for the variables, whether default or passed in, are stored in the output file (wigner.xsil).</p>
<div class="highlight-xpdeint"><div class="highlight"><pre><span class="nt"><info></span>
Script compiled with XMDS2 version 2.1 "Happy Mollusc" (r2680)
See http://www.xmds.org for more information.
Variables that can be specified on the command line:
Command line argument omega = 1.000000e-01
Command line argument alpha_0 = 2.500000e+00
Command line argument absorb = 8.000000e+00
Command line argument width = 3.000000e-01
Command line argument Uint_hbar = 0.000000e+00
<span class="nt"></info></span>
</pre></div>
</div>
<p id="index-12">Finally, note the shorthand used in the output group</p>
<div class="highlight-xpdeint"><div class="highlight"><pre><span class="cp"><![CDATA[</span>
<span class="n">_SAMPLE_COMPLEX</span><span class="p">(</span><span class="n">W</span><span class="p">);</span>
<span class="cp">]]></span>
</pre></div>
</div>
<p>which is short for</p>
<div class="highlight-xpdeint"><div class="highlight"><pre><span class="cp"><![CDATA[</span>
<span class="n">WR</span> <span class="o">=</span> <span class="n">W</span><span class="p">.</span><span class="n">Re</span><span class="p">();</span>
<span class="n">WI</span> <span class="o">=</span> <span class="n">W</span><span class="p">.</span><span class="n">Im</span><span class="p">();</span>
<span class="cp">]]></span>
</pre></div>
</div>
</div>
<div class="section" id="finding-the-ground-state-of-a-bec-continuous-renormalisation">
<span id="groundstatebec"></span><h2>Finding the Ground State of a BEC (continuous renormalisation)<a class="headerlink" href="#finding-the-ground-state-of-a-bec-continuous-renormalisation" title="Permalink to this headline">¶</a></h2>
<p>This simulation solves another partial differential equation, but introduces several powerful new features in XMDS2. The nominal problem is the calculation of the lowest energy eigenstate of a non-linear Schrödinger equation:</p>
<div class="math">
\[\frac{\partial \phi}{\partial t} = i \left[\frac{1}{2}\frac{\partial^2}{\partial y^2} - V(y) - U_{int}|\phi|^2\right]\phi\]</div>
<p>which can be found by evolving the above equation in imaginary time while keeping the normalisation constant. This causes eigenstates to exponentially decay at the rate of their eigenvalue, so after a short time only the state with the lowest eigenvalue remains. The evolution equation is straightforward:</p>
<div class="math">
\[\frac{\partial \phi}{\partial t} = \left[\frac{1}{2}\frac{\partial^2}{\partial y^2} - V(y) - U_{int}|\phi|^2\right]\phi\]</div>
<p>but we will need to use new XMDS2 features to manage the normalisation of the function <span class="math">\(\phi(y,t)\)</span>. The normalisation for a non-linear Schrödinger equation is given by <span class="math">\(\int dy |\phi(y,t)|^2 = N_{particles}\)</span>, where <span class="math">\(N_{particles}\)</span> is the number of particles described by the wavefunction.</p>
<p>The code for this simulation can be found in <tt class="docutils literal"><span class="pre">examples/groundstate_workedexamples.xmds</span></tt>:</p>
<div class="highlight-xpdeint"><div class="highlight"><pre><span class="nt"><simulation</span> <span class="na">xmds-version=</span><span class="s">"2"</span><span class="nt">></span>
<span class="nt"><name></span>groundstate<span class="nt"></name></span>
<span class="nt"><author></span>Joe Hope<span class="nt"></author></span>
<span class="nt"><description></span>
Calculate the ground state of the non-linear Schrodinger equation in a harmonic magnetic trap.
This is done by evolving it in imaginary time while re-normalising each timestep.
<span class="nt"></description></span>
<span class="nt"><features></span>
<span class="nt"><auto_vectorise</span> <span class="nt">/></span>
<span class="nt"><benchmark</span> <span class="nt">/></span>
<span class="nt"><bing</span> <span class="nt">/></span>
<span class="nt"><fftw</span> <span class="na">plan=</span><span class="s">"exhaustive"</span> <span class="nt">/></span>
<span class="nt"><globals></span>
<span class="cp"><![CDATA[</span>
<span class="k">const</span> <span class="nf">real</span> <span class="n">Uint</span> <span class="o">=</span> <span class="mf">2.0</span><span class="p">;</span>
<span class="k">const</span> <span class="nf">real</span> <span class="n">Nparticles</span> <span class="o">=</span> <span class="mf">5.0</span><span class="p">;</span>
<span class="cp">]]></span>
<span class="nt"></globals></span>
<span class="nt"></features></span>
<span class="nt"><geometry></span>
<span class="nt"><propagation_dimension></span> t <span class="nt"></propagation_dimension></span>
<span class="nt"><transverse_dimensions></span>
<span class="nt"><dimension</span> <span class="na">name=</span><span class="s">"y"</span> <span class="na">lattice=</span><span class="s">"256"</span> <span class="na">domain=</span><span class="s">"(-15.0, 15.0)"</span> <span class="nt">/></span>
<span class="nt"></transverse_dimensions></span>
<span class="nt"></geometry></span>
<span class="nt"><vector</span> <span class="na">name=</span><span class="s">"potential"</span> <span class="na">initial_basis=</span><span class="s">"y"</span> <span class="na">type=</span><span class="s">"real"</span><span class="nt">></span>
<span class="nt"><components></span> V1 <span class="nt"></components></span>
<span class="nt"><initialisation></span>
<span class="cp"><![CDATA[</span>
<span class="n">V1</span> <span class="o">=</span> <span class="mf">0.5</span><span class="o">*</span><span class="n">y</span><span class="o">*</span><span class="n">y</span><span class="p">;</span>
<span class="cp">]]></span>
<span class="nt"></initialisation></span>
<span class="nt"></vector></span>
<span class="nt"><vector</span> <span class="na">name=</span><span class="s">"wavefunction"</span> <span class="na">initial_basis=</span><span class="s">"y"</span> <span class="na">type=</span><span class="s">"complex"</span><span class="nt">></span>
<span class="nt"><components></span> phi <span class="nt"></components></span>
<span class="nt"><initialisation></span>
<span class="cp"><![CDATA[</span>
<span class="k">if</span> <span class="p">(</span><span class="n">fabs</span><span class="p">(</span><span class="n">y</span><span class="p">)</span> <span class="o"><</span> <span class="mf">3.0</span><span class="p">)</span> <span class="p">{</span>
<span class="n">phi</span> <span class="o">=</span> <span class="mf">1.0</span><span class="p">;</span>
<span class="c1">// This will be automatically normalised later</span>
<span class="p">}</span> <span class="k">else</span> <span class="p">{</span>
<span class="n">phi</span> <span class="o">=</span> <span class="mf">0.0</span><span class="p">;</span>
<span class="p">}</span>
<span class="cp">]]></span>
<span class="nt"></initialisation></span>
<span class="nt"></vector></span>
<span class="nt"><computed_vector</span> <span class="na">name=</span><span class="s">"normalisation"</span> <span class="na">dimensions=</span><span class="s">""</span> <span class="na">type=</span><span class="s">"real"</span><span class="nt">></span>
<span class="nt"><components></span> Ncalc <span class="nt"></components></span>
<span class="nt"><evaluation></span>
<span class="nt"><dependencies</span> <span class="na">basis=</span><span class="s">"y"</span><span class="nt">></span>wavefunction<span class="nt"></dependencies></span>
<span class="cp"><![CDATA[</span>
<span class="c1">// Calculate the current normalisation of the wave function.</span>
<span class="n">Ncalc</span> <span class="o">=</span> <span class="nf">mod2</span><span class="p">(</span><span class="n">phi</span><span class="p">);</span>
<span class="cp">]]></span>
<span class="nt"></evaluation></span>
<span class="nt"></computed_vector></span>
<span class="nt"><sequence></span>
<span class="nt"><filter></span>
<span class="cp"><![CDATA[</span>
<span class="n">printf</span><span class="p">(</span><span class="s">"Hello world from a filter segment!</span><span class="se">\n</span><span class="s">"</span><span class="p">);</span>
<span class="cp">]]></span>
<span class="nt"></filter></span>
<span class="nt"><filter></span>
<span class="nt"><dependencies></span>normalisation wavefunction<span class="nt"></dependencies></span>
<span class="cp"><![CDATA[</span>
<span class="n">phi</span> <span class="o">*=</span> <span class="n">sqrt</span><span class="p">(</span><span class="n">Nparticles</span><span class="o">/</span><span class="n">Ncalc</span><span class="p">);</span>
<span class="cp">]]></span>
<span class="nt"></filter></span>
<span class="nt"><integrate</span> <span class="na">algorithm=</span><span class="s">"ARK45"</span> <span class="na">interval=</span><span class="s">"1.0"</span> <span class="na">steps=</span><span class="s">"4000"</span> <span class="na">tolerance=</span><span class="s">"1e-10"</span><span class="nt">></span>
<span class="nt"><samples></span>25 4000<span class="nt"></samples></span>
<span class="nt"><filters</span> <span class="na">where=</span><span class="s">"step end"</span><span class="nt">></span>
<span class="nt"><filter></span>
<span class="nt"><dependencies></span>wavefunction normalisation<span class="nt"></dependencies></span>
<span class="cp"><![CDATA[</span>
<span class="c1">// Correct normalisation of the wavefunction</span>
<span class="n">phi</span> <span class="o">*=</span> <span class="n">sqrt</span><span class="p">(</span><span class="n">Nparticles</span><span class="o">/</span><span class="n">Ncalc</span><span class="p">);</span>
<span class="cp">]]></span>
<span class="nt"></filter></span>
<span class="nt"></filters></span>
<span class="nt"><operators></span>
<span class="nt"><operator</span> <span class="na">kind=</span><span class="s">"ip"</span><span class="nt">></span>
<span class="nt"><operator_names></span>T<span class="nt"></operator_names></span>
<span class="cp"><![CDATA[</span>
<span class="n">T</span> <span class="o">=</span> <span class="o">-</span><span class="mf">0.5</span><span class="o">*</span><span class="n">ky</span><span class="o">*</span><span class="n">ky</span><span class="p">;</span>
<span class="cp">]]></span>
<span class="nt"></operator></span>
<span class="nt"><integration_vectors></span>wavefunction<span class="nt"></integration_vectors></span>
<span class="nt"><dependencies></span>potential<span class="nt"></dependencies></span>
<span class="cp"><![CDATA[</span>
<span class="n">dphi_dt</span> <span class="o">=</span> <span class="n">T</span><span class="p">[</span><span class="n">phi</span><span class="p">]</span> <span class="o">-</span> <span class="p">(</span><span class="n">V1</span> <span class="o">+</span> <span class="n">Uint</span><span class="o">*</span><span class="nf">mod2</span><span class="p">(</span><span class="n">phi</span><span class="p">))</span><span class="o">*</span><span class="n">phi</span><span class="p">;</span>
<span class="cp">]]></span>
<span class="nt"></operators></span>
<span class="nt"></integrate></span>
<span class="nt"><breakpoint</span> <span class="na">filename=</span><span class="s">"groundstate_break.xsil"</span><span class="nt">></span>
<span class="nt"><dependencies</span> <span class="na">basis=</span><span class="s">"ky"</span><span class="nt">></span>wavefunction <span class="nt"></dependencies></span>
<span class="nt"></breakpoint></span>
<span class="nt"></sequence></span>
<span class="nt"><output></span>
<span class="nt"><sampling_group</span> <span class="na">basis=</span><span class="s">"y"</span> <span class="na">initial_sample=</span><span class="s">"yes"</span><span class="nt">></span>
<span class="nt"><moments></span>norm_dens<span class="nt"></moments></span>
<span class="nt"><dependencies></span>wavefunction normalisation<span class="nt"></dependencies></span>
<span class="cp"><![CDATA[</span>
<span class="n">norm_dens</span> <span class="o">=</span> <span class="nf">mod2</span><span class="p">(</span><span class="n">phi</span><span class="p">);</span>
<span class="cp">]]></span>
<span class="nt"></sampling_group></span>
<span class="nt"><sampling_group</span> <span class="na">initial_sample=</span><span class="s">"yes"</span><span class="nt">></span>
<span class="nt"><moments></span>norm<span class="nt"></moments></span>
<span class="nt"><dependencies></span>normalisation<span class="nt"></dependencies></span>
<span class="cp"><![CDATA[</span>
<span class="n">norm</span> <span class="o">=</span> <span class="n">Ncalc</span><span class="p">;</span>
<span class="cp">]]></span>
<span class="nt"></sampling_group></span>
<span class="nt"></output></span>
<span class="nt"></simulation></span>
</pre></div>
</div>
<p>We have used the <tt class="docutils literal"><span class="pre">plan="exhaustive"</span></tt> option in the <tt class="docutils literal"><span class="pre"><fftw></span></tt> element to ensure that the absolute fastest transform method is found. Because the FFTW package stores the results of its tests (by default in the ~/.xmds/wisdom directory), this option does not cause significant computational overhead, except perhaps on the very first run of a new program.</p>
<p id="index-13">This simulation introduces the first example of a very powerful feature in XMDS2: the <tt class="docutils literal"><span class="pre"><computed_vector></span></tt> element. This has syntax like any other vector, including possible dependencies on other vectors, and an ability to be used in any element that can use vectors. The difference is that, much like noise vectors, computed vectors are recalculated each time they are required. This means that a computed vector can never be used as an integration vector, as its values are not stored. However, computed vectors allow a simple and efficient method of describing complicated functions of other vectors. Computed vectors may depend on other computed vectors, allowing for spectral filtering and other advanced options. See for example, the <a class="reference internal" href="advanced_topics.html#advancedtopics"><em>Advanced Topics</em></a> section on <a class="reference internal" href="advanced_topics.html#convolutions"><em>Convolutions and Fourier transforms</em></a>.</p>
<p>The difference between a computed vector and a stored vector is emphasised by the replacement of the <tt class="docutils literal"><span class="pre"><initialisation></span></tt> element with an <tt class="docutils literal"><span class="pre"><evaluation></span></tt> element. Apart from the name, they have virtually identical purpose and syntax.</p>
<div class="highlight-xpdeint"><div class="highlight"><pre><span class="nt"><computed_vector</span> <span class="na">name=</span><span class="s">"normalisation"</span> <span class="na">dimensions=</span><span class="s">""</span> <span class="na">type=</span><span class="s">"real"</span><span class="nt">></span>
<span class="nt"><components></span> Ncalc <span class="nt"></components></span>
<span class="nt"><evaluation></span>
<span class="nt"><dependencies</span> <span class="na">basis=</span><span class="s">"y"</span><span class="nt">></span>wavefunction<span class="nt"></dependencies></span>
<span class="cp"><![CDATA[</span>
<span class="c1">// Calculate the current normalisation of the wave function.</span>
<span class="n">Ncalc</span> <span class="o">=</span> <span class="nf">mod2</span><span class="p">(</span><span class="n">phi</span><span class="p">);</span>
<span class="cp">]]></span>
<span class="nt"></evaluation></span>
<span class="nt"></computed_vector></span>
</pre></div>
</div>
<p>Here, our computed vector has no transverse dimensions and depends on the components of “wavefunction”, so the extra transverse dimensions are integrated out. This code therefore integrates the square modulus of the field, and returns it in the variable “Ncalc”. This will be used below to renormalise the “phi” field. Before we examine that process, we have to introduce the <tt class="docutils literal"><span class="pre"><filter></span></tt> element.</p>
<p>The <tt class="docutils literal"><span class="pre"><filter></span></tt> element can be placed in the <tt class="docutils literal"><span class="pre"><sequence></span></tt> element, or inside <tt class="docutils literal"><span class="pre"><integrate></span></tt> elements as we will see next. Elements placed in the <tt class="docutils literal"><span class="pre"><sequence></span></tt> element are executed in the order they are found in the .xmds file. Filter elements place the included CDATA block directly into the generated program at the designated position. If the element does not contain any dependencies, like in our first example, then the code is placed alone:</p>
<div class="highlight-xpdeint"><div class="highlight"><pre><span class="nt"><filter></span>
<span class="cp"><![CDATA[</span>
<span class="n">printf</span><span class="p">(</span><span class="s">"Hello world from a filter segment!</span><span class="se">\n</span><span class="s">"</span><span class="p">);</span>
<span class="cp">]]></span>
<span class="nt"></filter></span>
</pre></div>
</div>
<p>This filter block merely prints a string into the output when the generated program is run. If the <tt class="docutils literal"><span class="pre"><filter></span></tt> element contains dependencies, then the variables defined in those vectors (or computed vectors, or noise vectors) will be available, and the CDATA block will be placed inside loops that run over all the transverse dimensions used by the included vectors. The second filter block in this example depends on both the “wavefunction” and “normalisation” vectors:</p>
<div class="highlight-xpdeint"><div class="highlight"><pre><span class="nt"><filter></span>
<span class="nt"><dependencies></span>normalisation wavefunction<span class="nt"></dependencies></span>
<span class="cp"><![CDATA[</span>
<span class="n">phi</span> <span class="o">*=</span> <span class="n">sqrt</span><span class="p">(</span><span class="n">Nparticles</span><span class="o">/</span><span class="n">Ncalc</span><span class="p">);</span>
<span class="cp">]]></span>
<span class="nt"></filter></span>
</pre></div>
</div>
<p>Since this filter depends on a vector with the transverse dimension “y”, this filter will execute for each point in “y”. This code multiplies the value of the field “phi” by the factor required to produce a normalised function in the sense that <span class="math">\(\int dy |\phi(y,t)|^2 = N_{particles}\)</span>.</p>
<p>The next usage of a <tt class="docutils literal"><span class="pre"><filter></span></tt> element in this program is inside the <tt class="docutils literal"><span class="pre"><integrate></span></tt> element, where all filters are placed inside a <tt class="docutils literal"><span class="pre"><filters></span></tt> element.</p>
<div class="highlight-xpdeint"><div class="highlight"><pre><span class="nt"><filters</span> <span class="na">where=</span><span class="s">"step end"</span><span class="nt">></span>
<span class="nt"><filter></span>
<span class="nt"><dependencies></span>wavefunction normalisation<span class="nt"></dependencies></span>
<span class="cp"><![CDATA[</span>
<span class="c1">// Correct normalisation of the wavefunction</span>
<span class="n">phi</span> <span class="o">*=</span> <span class="n">sqrt</span><span class="p">(</span><span class="n">Nparticles</span><span class="o">/</span><span class="n">Ncalc</span><span class="p">);</span>
<span class="cp">]]></span>
<span class="nt"></filter></span>
<span class="nt"></filters></span>
</pre></div>
</div>
<p>Filters placed in an integration block are applied each integration step. The “where” flag is used to determine whether the filter should be applied directly before or directly after each integration step. The default value for the where flag is <tt class="docutils literal"><span class="pre">where="step</span> <span class="pre">start"</span></tt>, but in this case we chose “step end” to make sure that the final output was normalised after the last integration step.</p>
<p>At the end of the sequence element we introduce the <tt class="docutils literal"><span class="pre"><breakpoint></span></tt> element. This serves two purposes. The first is a simple matter of convenience. Often when we manage our input and output from a simulation, we are interested solely in storing the exact state of our integration vectors. A breakpoint element does exactly that, storing the components of any vectors contained within, taking all the normal options of the <tt class="docutils literal"><span class="pre"><output></span></tt> element but not requiring any <tt class="docutils literal"><span class="pre"><sampling_group></span></tt> elements as that information is assumed.</p>
<div class="highlight-xpdeint"><div class="highlight"><pre><span class="nt"><breakpoint</span> <span class="na">filename=</span><span class="s">"groundstate_break.xsil"</span><span class="nt">></span>
<span class="nt"><dependencies</span> <span class="na">basis=</span><span class="s">"ky"</span><span class="nt">></span>wavefunction<span class="nt"></dependencies></span>
<span class="nt"></breakpoint></span>
</pre></div>
</div>
<p>If the filename argument is omitted, the output filenames are numbered sequentially. Any given <tt class="docutils literal"><span class="pre"><breakpoint></span></tt> element must only depend on vectors with identical dimensions.</p>
<p>This program begins with a very crude guess to the ground state, but it rapidly converges to the lowest eigenstate.</p>
<div class="figure align-center">
<img alt="_images/groundstateU2.png" src="_images/groundstateU2.png" />
<p class="caption">The shape of the ground state rapidly approaches the lowest eigenstate. For weak nonlinearities, it is nearly Gaussian.</p>
</div>
<div class="figure align-center">
<img alt="_images/groundstateU20.png" src="_images/groundstateU20.png" />
<p class="caption">When the nonlinear term is larger (<span class="math">\(U=20\)</span>), the ground state is wider and more parabolic.</p>
</div>
</div>
<div class="section" id="finding-the-ground-state-of-a-bec-again">
<span id="hermitegaussgroundstatebec"></span><h2>Finding the Ground State of a BEC again<a class="headerlink" href="#finding-the-ground-state-of-a-bec-again" title="Permalink to this headline">¶</a></h2>
<p>Here we repeat the same simulation as in the <a class="reference internal" href="#groundstatebec"><em>Finding the Ground State of a BEC (continuous renormalisation)</em></a> example, using a different transform basis. While spectral methods are very effective, and Fourier transforms are typically very efficient due to the Fast Fourier transform algorithm, it is often desirable to describe nonlocal evolution in bases other than the Fourier basis. The previous calculation was the Schrödinger equation with a harmonic potential and a nonlinear term. The eigenstates of such a system are known analytically to be Gaussians multiplied by the Hermite polynomials.</p>
<div class="math">
\[\left[-\frac{\hbar}{2 m}\frac{\partial^2}{\partial x^2} + \frac{1}{2}\omega^2 x^2\right]\phi_n(x) = E_n \phi_n(x)\]</div>
<p>where</p>
<div class="math">
\[\phi_n(x,t) = \sqrt{\frac{1}{2^n n!}} \left(\frac{m \omega}{\hbar \pi}\right)^\frac{1}{4} e^{-\frac{m \omega x^2}{2\hbar}} H_n\left(\sqrt{\frac{m \omega}{\hbar}x}\right),\;\;\;\;\;\;E_n = \left(n+\frac{1}{2}\right) \omega\]</div>
<p>where <span class="math">\(H_n(u)\)</span> are the physicist’s version of the Hermite polynomials. Rather than describing the derivatives as diagonal terms in Fourier space, we therefore have the option of describing the entire <span class="math">\(-\frac{\hbar}{2 m}\frac{\partial^2}{\partial x^2} + \frac{1}{2}\omega^2 x^2\)</span> term as a diagonal term in the Hermite-Gaussian basis. Here is an XMDS2 simulation that performs the integration in this basis. The following is a simplified version of the <tt class="docutils literal"><span class="pre">examples/hermitegauss_groundstate.xmds</span></tt> script.</p>
<div class="highlight-xpdeint"><div class="highlight"><pre><span class="nt"><simulation</span> <span class="na">xmds-version=</span><span class="s">"2"</span><span class="nt">></span>
<span class="nt"><name></span>hermitegauss_groundstate<span class="nt"></name></span>
<span class="nt"><author></span>Graham Dennis<span class="nt"></author></span>
<span class="nt"><description></span>
Solve for the groundstate of the Gross-Pitaevskii equation using the Hermite-Gauss basis.
<span class="nt"></description></span>
<span class="nt"><features></span>
<span class="nt"><benchmark</span> <span class="nt">/></span>
<span class="nt"><bing</span> <span class="nt">/></span>
<span class="nt"><validation</span> <span class="na">kind=</span><span class="s">"run-time"</span> <span class="nt">/></span>
<span class="nt"><globals></span>
<span class="cp"><![CDATA[</span>
<span class="k">const</span> <span class="nf">real</span> <span class="n">omegaz</span> <span class="o">=</span> <span class="mi">2</span><span class="o">*</span><span class="n">M_PI</span><span class="o">*</span><span class="mi">20</span><span class="p">;</span>
<span class="k">const</span> <span class="nf">real</span> <span class="n">omegarho</span> <span class="o">=</span> <span class="mi">2</span><span class="o">*</span><span class="n">M_PI</span><span class="o">*</span><span class="mi">200</span><span class="p">;</span>
<span class="k">const</span> <span class="nf">real</span> <span class="n">hbar</span> <span class="o">=</span> <span class="mf">1.05457148e-34</span><span class="p">;</span>
<span class="k">const</span> <span class="nf">real</span> <span class="n">M</span> <span class="o">=</span> <span class="mf">1.409539200000000e-25</span><span class="p">;</span>
<span class="k">const</span> <span class="nf">real</span> <span class="n">g</span> <span class="o">=</span> <span class="mf">9.8</span><span class="p">;</span>
<span class="k">const</span> <span class="nf">real</span> <span class="n">scatteringLength</span> <span class="o">=</span> <span class="mf">5.57e-9</span><span class="p">;</span>
<span class="k">const</span> <span class="nf">real</span> <span class="n">transverseLength</span> <span class="o">=</span> <span class="mf">1e-5</span><span class="p">;</span>
<span class="k">const</span> <span class="nf">real</span> <span class="n">Uint</span> <span class="o">=</span> <span class="mf">4.0</span><span class="o">*</span><span class="n">M_PI</span><span class="o">*</span><span class="n">hbar</span><span class="o">*</span><span class="n">hbar</span><span class="o">*</span><span class="n">scatteringLength</span><span class="o">/</span><span class="n">M</span><span class="o">/</span><span class="n">transverseLength</span><span class="o">/</span><span class="n">transverseLength</span><span class="p">;</span>
<span class="k">const</span> <span class="nf">real</span> <span class="n">Nparticles</span> <span class="o">=</span> <span class="mf">5.0e5</span><span class="p">;</span>
<span class="cm">/* offset constants */</span>
<span class="k">const</span> <span class="nf">real</span> <span class="n">EnergyOffset</span> <span class="o">=</span> <span class="mf">0.3</span><span class="o">*</span><span class="n">pow</span><span class="p">(</span><span class="n">pow</span><span class="p">(</span><span class="mf">3.0</span><span class="o">*</span><span class="n">Nparticles</span><span class="o">/</span><span class="mi">4</span><span class="o">*</span><span class="n">omegarho</span><span class="o">*</span><span class="n">Uint</span><span class="p">,</span><span class="mf">2.0</span><span class="p">)</span><span class="o">*</span><span class="n">M</span><span class="o">/</span><span class="mf">2.0</span><span class="p">,</span><span class="mi">1</span><span class="o">/</span><span class="mf">3.0</span><span class="p">);</span> <span class="c1">// 1D</span>
<span class="cp">]]></span>
<span class="nt"></globals></span>
<span class="nt"></features></span>
<span class="nt"><geometry></span>
<span class="nt"><propagation_dimension></span> t <span class="nt"></propagation_dimension></span>
<span class="nt"><transverse_dimensions></span>
<span class="nt"><dimension</span> <span class="na">name=</span><span class="s">"x"</span> <span class="na">lattice=</span><span class="s">"100"</span> <span class="na">length_scale=</span><span class="s">"sqrt(hbar/(M*omegarho))"</span> <span class="na">transform=</span><span class="s">"hermite-gauss"</span> <span class="nt">/></span>
<span class="nt"></transverse_dimensions></span>
<span class="nt"></geometry></span>
<span class="nt"><vector</span> <span class="na">name=</span><span class="s">"wavefunction"</span> <span class="na">initial_basis=</span><span class="s">"x"</span> <span class="na">type=</span><span class="s">"complex"</span><span class="nt">></span>
<span class="nt"><components></span> phi <span class="nt"></components></span>
<span class="nt"><initialisation></span>
<span class="cp"><![CDATA[</span>
<span class="n">phi</span> <span class="o">=</span> <span class="n">sqrt</span><span class="p">(</span><span class="n">Nparticles</span><span class="p">)</span> <span class="o">*</span> <span class="n">pow</span><span class="p">(</span><span class="n">M</span><span class="o">*</span><span class="n">omegarho</span><span class="o">/</span><span class="p">(</span><span class="n">hbar</span><span class="o">*</span><span class="n">M_PI</span><span class="p">),</span> <span class="mf">0.25</span><span class="p">)</span> <span class="o">*</span> <span class="n">exp</span><span class="p">(</span><span class="o">-</span><span class="mf">0.5</span><span class="o">*</span><span class="p">(</span><span class="n">M</span><span class="o">*</span><span class="n">omegarho</span><span class="o">/</span><span class="n">hbar</span><span class="p">)</span><span class="o">*</span><span class="n">x</span><span class="o">*</span><span class="n">x</span><span class="p">);</span>
<span class="cp">]]></span>
<span class="nt"></initialisation></span>
<span class="nt"></vector></span>
<span class="nt"><computed_vector</span> <span class="na">name=</span><span class="s">"normalisation"</span> <span class="na">dimensions=</span><span class="s">""</span> <span class="na">type=</span><span class="s">"real"</span><span class="nt">></span>
<span class="nt"><components></span> Ncalc <span class="nt"></components></span>
<span class="nt"><evaluation></span>
<span class="nt"><dependencies</span> <span class="na">basis=</span><span class="s">"x"</span><span class="nt">></span>wavefunction<span class="nt"></dependencies></span>
<span class="cp"><![CDATA[</span>
<span class="c1">// Calculate the current normalisation of the wave function.</span>
<span class="n">Ncalc</span> <span class="o">=</span> <span class="nf">mod2</span><span class="p">(</span><span class="n">phi</span><span class="p">);</span>
<span class="cp">]]></span>
<span class="nt"></evaluation></span>
<span class="nt"></computed_vector></span>
<span class="nt"><sequence></span>
<span class="nt"><integrate</span> <span class="na">algorithm=</span><span class="s">"ARK45"</span> <span class="na">interval=</span><span class="s">"1.0e-2"</span> <span class="na">steps=</span><span class="s">"4000"</span> <span class="na">tolerance=</span><span class="s">"1e-10"</span><span class="nt">></span>
<span class="nt"><samples></span>100 100<span class="nt"></samples></span>
<span class="nt"><filters></span>
<span class="nt"><filter></span>
<span class="nt"><dependencies></span>wavefunction normalisation<span class="nt"></dependencies></span>
<span class="cp"><![CDATA[</span>
<span class="c1">// Correct normalisation of the wavefunction</span>
<span class="n">phi</span> <span class="o">*=</span> <span class="n">sqrt</span><span class="p">(</span><span class="n">Nparticles</span><span class="o">/</span><span class="n">Ncalc</span><span class="p">);</span>
<span class="cp">]]></span>
<span class="nt"></filter></span>
<span class="nt"></filters></span>
<span class="nt"><operators></span>
<span class="nt"><operator</span> <span class="na">kind=</span><span class="s">"ip"</span> <span class="na">type=</span><span class="s">"real"</span><span class="nt">></span>
<span class="nt"><operator_names></span>L<span class="nt"></operator_names></span>
<span class="cp"><![CDATA[</span>
<span class="n">L</span> <span class="o">=</span> <span class="n">EnergyOffset</span><span class="o">/</span><span class="n">hbar</span> <span class="o">-</span> <span class="p">(</span><span class="n">nx</span> <span class="o">+</span> <span class="mf">0.5</span><span class="p">)</span><span class="o">*</span><span class="n">omegarho</span><span class="p">;</span>
<span class="cp">]]></span>
<span class="nt"></operator></span>
<span class="nt"><integration_vectors></span>wavefunction<span class="nt"></integration_vectors></span>
<span class="cp"><![CDATA[</span>
<span class="n">dphi_dt</span> <span class="o">=</span> <span class="n">L</span><span class="p">[</span><span class="n">phi</span><span class="p">]</span> <span class="o">-</span> <span class="n">Uint</span><span class="o">/</span><span class="n">hbar</span><span class="o">*</span><span class="nf">mod2</span><span class="p">(</span><span class="n">phi</span><span class="p">)</span><span class="o">*</span><span class="n">phi</span><span class="p">;</span>
<span class="cp">]]></span>
<span class="nt"></operators></span>
<span class="nt"></integrate></span>
<span class="nt"><filter></span>
<span class="nt"><dependencies></span>normalisation wavefunction<span class="nt"></dependencies></span>
<span class="cp"><![CDATA[</span>
<span class="n">phi</span> <span class="o">*=</span> <span class="n">sqrt</span><span class="p">(</span><span class="n">Nparticles</span><span class="o">/</span><span class="n">Ncalc</span><span class="p">);</span>
<span class="cp">]]></span>
<span class="nt"></filter></span>
<span class="nt"><breakpoint</span> <span class="na">filename=</span><span class="s">"hermitegauss_groundstate_break.xsil"</span> <span class="na">format=</span><span class="s">"ascii"</span><span class="nt">></span>
<span class="nt"><dependencies</span> <span class="na">basis=</span><span class="s">"nx"</span><span class="nt">></span>wavefunction<span class="nt"></dependencies></span>
<span class="nt"></breakpoint></span>
<span class="nt"></sequence></span>
<span class="nt"><output></span>
<span class="nt"><sampling_group</span> <span class="na">basis=</span><span class="s">"x"</span> <span class="na">initial_sample=</span><span class="s">"yes"</span><span class="nt">></span>
<span class="nt"><moments></span>dens<span class="nt"></moments></span>
<span class="nt"><dependencies></span>wavefunction<span class="nt"></dependencies></span>
<span class="cp"><![CDATA[</span>
<span class="n">dens</span> <span class="o">=</span> <span class="nf">mod2</span><span class="p">(</span><span class="n">phi</span><span class="p">);</span>
<span class="cp">]]></span>
<span class="nt"></sampling_group></span>
<span class="nt"><sampling_group</span> <span class="na">basis=</span><span class="s">"kx"</span> <span class="na">initial_sample=</span><span class="s">"yes"</span><span class="nt">></span>
<span class="nt"><moments></span>dens<span class="nt"></moments></span>
<span class="nt"><dependencies></span>wavefunction<span class="nt"></dependencies></span>
<span class="cp"><![CDATA[</span>
<span class="n">dens</span> <span class="o">=</span> <span class="nf">mod2</span><span class="p">(</span><span class="n">phi</span><span class="p">);</span>
<span class="cp">]]></span>
<span class="nt"></sampling_group></span>
<span class="nt"></output></span>
<span class="nt"></simulation></span>
</pre></div>
</div>
<p>The major difference in this simulation code, aside from the switch back from dimensionless units, is the new transverse dimension type in the <tt class="docutils literal"><span class="pre"><geometry></span></tt> element.</p>
<div class="highlight-xpdeint"><div class="highlight"><pre><span class="nt"><dimension</span> <span class="na">name=</span><span class="s">"x"</span> <span class="na">lattice=</span><span class="s">"100"</span> <span class="na">length_scale=</span><span class="s">"sqrt(hbar/(M*omegarho))"</span> <span class="na">transform=</span><span class="s">"hermite-gauss"</span> <span class="nt">/></span>
</pre></div>
</div>
<p>We have explicitly defined the “transform” option, which by default expects the Fourier transform. The <tt class="docutils literal"><span class="pre">transform="hermite-gauss"</span></tt> option requires the ‘mpmath’ package installed, just as Fourier transforms require the FFTW package to be installed. The “lattice” option details the number of hermite-Gaussian eigenstates to include, and automatically starts from the zeroth order polynomial and increases. The number of hermite-Gaussian modes fully determines the irregular spatial grid up to an overall scale given by the <tt class="docutils literal"><span class="pre">length_scale</span></tt> parameter.</p>
<p>The <tt class="docutils literal"><span class="pre">length_scale="sqrt(hbar/(M*omegarho))"</span></tt> option requires a real number, but since this script defines it in terms of variables, XMDS2 is unable to verify that the resulting function is real-valued at the time of generating the code. XMDS2 will therefore fail to compile this program without the feature:</p>
<div class="highlight-xpdeint"><div class="highlight"><pre><span class="nt"><validation</span> <span class="na">kind=</span><span class="s">"run-time"</span> <span class="nt">/></span>
</pre></div>
</div>
<p>which disables many of these checks at the time of writing the C-code.</p>
</div>
<div class="section" id="multi-component-schrodinger-equation">
<span id="dmultistatese"></span><h2>Multi-component Schrödinger equation<a class="headerlink" href="#multi-component-schrodinger-equation" title="Permalink to this headline">¶</a></h2>
<p>This example demonstrates a simple method for doing matrix calculations in XMDS2. We are solving the multi-component PDE</p>
<div class="math">
\[\frac{\partial \phi_j(x,y)}{\partial t} = \frac{i}{2}\left(\frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial y^2}\right)\phi_j(x,y) - i U(x,y) \sum_k V_{j k}\phi_k(x,y)\]</div>
<p>where the last term is more commonly written as a matrix multiplication. Writing this term out explicitly is feasible for a small number of components, but when the number of components becomes large, or perhaps <span class="math">\(V_{j k}\)</span> should be precomputed for efficiency reasons, it is useful to be able to perform this sum over the integer dimensions automatically. This example show how this can be done naturally using a computed vector. The XMDS2 script is as follows:</p>
<div class="highlight-xpdeint"><div class="highlight"><pre><span class="nt"><simulation</span> <span class="na">xmds-version=</span><span class="s">"2"</span><span class="nt">></span>
<span class="nt"><name></span>2DMSse<span class="nt"></name></span>
<span class="nt"><author></span>Joe Hope<span class="nt"></author></span>
<span class="nt"><description></span>
Schroedinger equation for multiple internal states in two spatial dimensions.
<span class="nt"></description></span>
<span class="nt"><features></span>
<span class="nt"><benchmark</span> <span class="nt">/></span>
<span class="nt"><bing</span> <span class="nt">/></span>
<span class="nt"><fftw</span> <span class="na">plan=</span><span class="s">"patient"</span> <span class="nt">/></span>
<span class="nt"><openmp</span> <span class="nt">/></span>
<span class="nt"><auto_vectorise</span> <span class="nt">/></span>
<span class="nt"></features></span>
<span class="nt"><geometry></span>
<span class="nt"><propagation_dimension></span> t <span class="nt"></propagation_dimension></span>
<span class="nt"><transverse_dimensions></span>
<span class="nt"><dimension</span> <span class="na">name=</span><span class="s">"x"</span> <span class="na">lattice=</span><span class="s">"32"</span> <span class="na">domain=</span><span class="s">"(-6, 6)"</span> <span class="nt">/></span>
<span class="nt"><dimension</span> <span class="na">name=</span><span class="s">"y"</span> <span class="na">lattice=</span><span class="s">"32"</span> <span class="na">domain=</span><span class="s">"(-6, 6)"</span> <span class="nt">/></span>
<span class="nt"><dimension</span> <span class="na">name=</span><span class="s">"j"</span> <span class="na">type=</span><span class="s">"integer"</span> <span class="na">lattice=</span><span class="s">"2"</span> <span class="na">domain=</span><span class="s">"(0,1)"</span> <span class="na">aliases=</span><span class="s">"k"</span><span class="nt">/></span>
<span class="nt"></transverse_dimensions></span>
<span class="nt"></geometry></span>
<span class="nt"><vector</span> <span class="na">name=</span><span class="s">"wavefunction"</span> <span class="na">type=</span><span class="s">"complex"</span> <span class="na">dimensions=</span><span class="s">"x y j"</span><span class="nt">></span>
<span class="nt"><components></span> phi <span class="nt"></components></span>
<span class="nt"><initialisation></span>
<span class="cp"><![CDATA[</span>
<span class="n">phi</span> <span class="o">=</span> <span class="n">j</span><span class="o">*</span><span class="n">sqrt</span><span class="p">(</span><span class="mi">2</span><span class="o">/</span><span class="n">sqrt</span><span class="p">(</span><span class="n">M_PI</span><span class="o">/</span><span class="mi">2</span><span class="p">))</span><span class="o">*</span><span class="n">exp</span><span class="p">(</span><span class="o">-</span><span class="p">(</span><span class="n">x</span><span class="o">*</span><span class="n">x</span><span class="o">+</span><span class="n">y</span><span class="o">*</span><span class="n">y</span><span class="p">)</span><span class="o">/</span><span class="mi">4</span><span class="p">)</span><span class="o">*</span><span class="n">exp</span><span class="p">(</span><span class="kc">i</span><span class="o">*</span><span class="mf">0.1</span><span class="o">*</span><span class="n">x</span><span class="p">);</span>
<span class="cp">]]></span>
<span class="nt"></initialisation></span>
<span class="nt"></vector></span>
<span class="nt"><vector</span> <span class="na">name=</span><span class="s">"spatialInteraction"</span> <span class="na">type=</span><span class="s">"real"</span> <span class="na">dimensions=</span><span class="s">"x y"</span><span class="nt">></span>
<span class="nt"><components></span> U <span class="nt"></components></span>
<span class="nt"><initialisation></span>
<span class="cp"><![CDATA[</span>
<span class="n">U</span><span class="o">=</span><span class="n">exp</span><span class="p">(</span><span class="o">-</span><span class="p">(</span><span class="n">x</span><span class="o">*</span><span class="n">x</span><span class="o">+</span><span class="n">y</span><span class="o">*</span><span class="n">y</span><span class="p">)</span><span class="o">/</span><span class="mi">4</span><span class="p">);</span>
<span class="cp">]]></span>
<span class="nt"></initialisation></span>
<span class="nt"></vector></span>
<span class="nt"><vector</span> <span class="na">name=</span><span class="s">"internalInteraction"</span> <span class="na">type=</span><span class="s">"real"</span> <span class="na">dimensions=</span><span class="s">"j k"</span><span class="nt">></span>
<span class="nt"><components></span> V <span class="nt"></components></span>
<span class="nt"><initialisation></span>
<span class="cp"><![CDATA[</span>
<span class="n">V</span><span class="o">=</span><span class="mi">3</span><span class="o">*</span><span class="p">(</span><span class="n">j</span><span class="o">*</span><span class="p">(</span><span class="mi">1</span><span class="o">-</span><span class="n">k</span><span class="p">)</span><span class="o">+</span><span class="p">(</span><span class="mi">1</span><span class="o">-</span><span class="n">j</span><span class="p">)</span><span class="o">*</span><span class="n">k</span><span class="p">);</span>
<span class="cp">]]></span>
<span class="nt"></initialisation></span>
<span class="nt"></vector></span>
<span class="nt"><computed_vector</span> <span class="na">name=</span><span class="s">"coupling"</span> <span class="na">dimensions=</span><span class="s">"x y j"</span> <span class="na">type=</span><span class="s">"complex"</span><span class="nt">></span>
<span class="nt"><components></span>
VPhi
<span class="nt"></components></span>
<span class="nt"><evaluation></span>
<span class="nt"><dependencies</span> <span class="na">basis=</span><span class="s">"x y j k"</span><span class="nt">></span>internalInteraction wavefunction<span class="nt"></dependencies></span>
<span class="cp"><![CDATA[</span>
<span class="c1">// Calculate the current normalisation of the wave function.</span>
<span class="n">VPhi</span> <span class="o">=</span> <span class="n">V</span><span class="o">*</span><span class="n">phi</span><span class="p">(</span><span class="n">j</span> <span class="o">=></span> <span class="n">k</span><span class="p">);</span>
<span class="cp">]]></span>
<span class="nt"></evaluation></span>
<span class="nt"></computed_vector></span>
<span class="nt"><sequence></span>
<span class="nt"><integrate</span> <span class="na">algorithm=</span><span class="s">"ARK45"</span> <span class="na">interval=</span><span class="s">"2.0"</span> <span class="na">tolerance=</span><span class="s">"1e-7"</span><span class="nt">></span>
<span class="nt"><samples></span>20 100<span class="nt"></samples></span>
<span class="nt"><operators></span>
<span class="nt"><integration_vectors></span>wavefunction<span class="nt"></integration_vectors></span>
<span class="nt"><operator</span> <span class="na">kind=</span><span class="s">"ip"</span> <span class="na">dimensions=</span><span class="s">"x"</span><span class="nt">></span>
<span class="nt"><operator_names></span>Lx<span class="nt"></operator_names></span>
<span class="cp"><![CDATA[</span>
<span class="n">Lx</span> <span class="o">=</span> <span class="o">-</span><span class="kc">i</span><span class="o">*</span><span class="n">kx</span><span class="o">*</span><span class="n">kx</span><span class="o">*</span><span class="mf">0.5</span><span class="p">;</span>
<span class="cp">]]></span>
<span class="nt"></operator></span>
<span class="nt"><operator</span> <span class="na">kind=</span><span class="s">"ip"</span> <span class="na">dimensions=</span><span class="s">"y"</span><span class="nt">></span>
<span class="nt"><operator_names></span>Ly<span class="nt"></operator_names></span>
<span class="cp"><![CDATA[</span>
<span class="n">Ly</span> <span class="o">=</span> <span class="o">-</span><span class="kc">i</span><span class="o">*</span><span class="n">ky</span><span class="o">*</span><span class="n">ky</span><span class="o">*</span><span class="mf">0.5</span><span class="p">;</span>
<span class="cp">]]></span>
<span class="nt"></operator></span>
<span class="cp"><![CDATA[</span>
<span class="n">dphi_dt</span> <span class="o">=</span> <span class="n">Lx</span><span class="p">[</span><span class="n">phi</span><span class="p">]</span> <span class="o">+</span> <span class="n">Ly</span><span class="p">[</span><span class="n">phi</span><span class="p">]</span> <span class="o">-</span><span class="kc">i</span><span class="o">*</span><span class="n">U</span><span class="o">*</span><span class="n">VPhi</span><span class="p">;</span>
<span class="cp">]]></span>
<span class="nt"><dependencies></span>spatialInteraction coupling<span class="nt"></dependencies></span>
<span class="nt"></operators></span>
<span class="nt"></integrate></span>
<span class="nt"></sequence></span>
<span class="nt"><output></span>
<span class="nt"><sampling_group</span> <span class="na">basis=</span><span class="s">"x y j"</span> <span class="na">initial_sample=</span><span class="s">"yes"</span><span class="nt">></span>
<span class="nt"><moments></span>density<span class="nt"></moments></span>
<span class="nt"><dependencies></span>wavefunction<span class="nt"></dependencies></span>
<span class="cp"><![CDATA[</span>
<span class="n">density</span> <span class="o">=</span> <span class="nf">mod2</span><span class="p">(</span><span class="n">phi</span><span class="p">);</span>
<span class="cp">]]></span>
<span class="nt"></sampling_group></span>
<span class="nt"><sampling_group</span> <span class="na">basis=</span><span class="s">"x(0) y(0) j"</span> <span class="na">initial_sample=</span><span class="s">"yes"</span><span class="nt">></span>
<span class="nt"><moments></span>normalisation<span class="nt"></moments></span>
<span class="nt"><dependencies></span>wavefunction<span class="nt"></dependencies></span>
<span class="cp"><![CDATA[</span>
<span class="n">normalisation</span> <span class="o">=</span> <span class="nf">mod2</span><span class="p">(</span><span class="n">phi</span><span class="p">);</span>
<span class="cp">]]></span>
<span class="nt"></sampling_group></span>
<span class="nt"></output></span>
<span class="nt"></simulation></span>
</pre></div>
</div>
<p id="index-14">The only truly new feature in this script is the “aliases” option on a dimension. The integer-valued dimension in this script indexes the components of the PDE (in this case only two). The <span class="math">\(V_{j k}\)</span> term is required to be a square array of dimension of this number of components. If we wrote the k-index of <span class="math">\(V_{j k}\)</span> using a separate <tt class="docutils literal"><span class="pre"><dimension></span></tt> element, then we would not be enforcing the requirement that the matrix be square. Instead, we note that we will be using multiple ‘copies’ of the j-dimension by using the “aliases” tag.</p>
<div class="highlight-xpdeint"><div class="highlight"><pre><span class="nt"><dimension</span> <span class="na">name=</span><span class="s">"j"</span> <span class="na">type=</span><span class="s">"integer"</span> <span class="na">lattice=</span><span class="s">"2"</span> <span class="na">domain=</span><span class="s">"(0,1)"</span> <span class="na">aliases=</span><span class="s">"k"</span><span class="nt">/></span>
</pre></div>
</div>
<p>This means that we can use the index “k”, which will have exactly the same properties as the “j” index. This is used to define the “V” function in the “internalInteraction” vector. Now, just as we use a computed vector to perform an integration over our fields, we use a computed vector to calculate the sum.</p>
<div class="highlight-xpdeint"><div class="highlight"><pre><span class="nt"><computed_vector</span> <span class="na">name=</span><span class="s">"coupling"</span> <span class="na">dimensions=</span><span class="s">"x y j"</span> <span class="na">type=</span><span class="s">"complex"</span><span class="nt">></span>
<span class="nt"><components></span>
VPhi
<span class="nt"></components></span>
<span class="nt"><evaluation></span>
<span class="nt"><dependencies</span> <span class="na">basis=</span><span class="s">"x y j k"</span><span class="nt">></span>internalInteraction wavefunction<span class="nt"></dependencies></span>
<span class="cp"><![CDATA[</span>
<span class="c1">// Calculate the current normalisation of the wave function.</span>
<span class="n">VPhi</span> <span class="o">=</span> <span class="n">V</span><span class="o">*</span><span class="n">phi</span><span class="p">(</span><span class="n">j</span> <span class="o">=></span> <span class="n">k</span><span class="p">);</span>
<span class="cp">]]></span>
<span class="nt"></evaluation></span>
<span class="nt"></computed_vector></span>
</pre></div>
</div>
<p>Since the output dimensions of the computed vector do not include a “k” index, this index is integrated. The volume element for this summation is the spacing between neighbouring values of “j”, and since this spacing is one, this integration is just a sum over k, as required.</p>
<p>This example also demonstrates an optimisation for the IP operators by separating the <span class="math">\(x\)</span> and <span class="math">\(y\)</span> parts of the operator (see <a class="reference internal" href="optimisation_hints.html#optimisingipoperators"><em>Optimising with the Interaction Picture (IP) operator</em></a>). This gives an approximately 30% speed improvement over the more straightforward implementation:</p>
<div class="highlight-xpdeint"><div class="highlight"><pre><span class="nt"><operator</span> <span class="na">kind=</span><span class="s">"ip"</span><span class="nt">></span>
<span class="nt"><operator_names></span>L<span class="nt"></operator_names></span>
<span class="cp"><![CDATA[</span>
<span class="n">L</span> <span class="o">=</span> <span class="o">-</span><span class="kc">i</span><span class="o">*</span><span class="p">(</span><span class="n">kx</span><span class="o">*</span><span class="n">kx</span> <span class="o">+</span> <span class="n">ky</span><span class="o">*</span><span class="n">ky</span><span class="p">)</span><span class="o">*</span><span class="mf">0.5</span><span class="p">;</span>
<span class="cp">]]></span>
<span class="nt"></operator></span>
<span class="cp"><![CDATA[</span>
<span class="n">dphi_dt</span> <span class="o">=</span> <span class="n">L</span><span class="p">[</span><span class="n">phi</span><span class="p">]</span> <span class="o">-</span> <span class="kc">i</span><span class="o">*</span><span class="n">U</span><span class="o">*</span><span class="n">VPhi</span><span class="p">;</span>
<span class="cp">]]></span>
</pre></div>
</div>
<p>By this point, we have introduced most of the important features in XMDS2. More details on other transform options and rarely used features can be found in the <a class="reference internal" href="advanced_topics.html#advancedtopics"><em>Advanced Topics</em></a> section.</p>
</div>
</div>
</div>
</div>
</div>
<div class="sphinxsidebar">
<div class="sphinxsidebarwrapper">
<p class="logo"><a href="index.html">
<img class="logo" src="_static/xmds_logo.png" alt="Logo"/>
</a></p>
<h3><a href="index.html">Table Of Contents</a></h3>
<ul>
<li><a class="reference internal" href="#">Worked Examples</a><ul>
<li><a class="reference internal" href="#the-nonlinear-schrodinger-equation">The nonlinear Schrödinger equation</a></li>
<li><a class="reference internal" href="#kubo-oscillator">Kubo Oscillator</a></li>
<li><a class="reference internal" href="#fibre-noise">Fibre Noise</a></li>
<li><a class="reference internal" href="#integer-dimensions">Integer Dimensions</a></li>
<li><a class="reference internal" href="#wigner-function">Wigner Function</a></li>
<li><a class="reference internal" href="#finding-the-ground-state-of-a-bec-continuous-renormalisation">Finding the Ground State of a BEC (continuous renormalisation)</a></li>
<li><a class="reference internal" href="#finding-the-ground-state-of-a-bec-again">Finding the Ground State of a BEC again</a></li>
<li><a class="reference internal" href="#multi-component-schrodinger-equation">Multi-component Schrödinger equation</a></li>
</ul>
</li>
</ul>
<div id="searchbox" style="display: none">
<h3>Quick search</h3>
<form class="search" action="search.html" method="get">
<input type="text" name="q" />
<input type="submit" value="Go" />
<input type="hidden" name="check_keywords" value="yes" />
<input type="hidden" name="area" value="default" />
</form>
<p class="searchtip" style="font-size: 90%">
Enter search terms or a module, class or function name.
</p>
</div>
<script type="text/javascript">$('#searchbox').show(0);</script>
</div>
</div>
<div class="clearer"></div>
</div>
<div class="related">
<h3>Navigation</h3>
<ul>
<li class="right" style="margin-right: 10px">
<a href="genindex.html" title="General Index"
>index</a></li>
<li class="right" >
<a href="reference_index.html" title="Reference section"
>next</a> |</li>
<li class="right" >
<a href="tutorial.html" title="Quickstart Tutorial"
>previous</a> |</li>
<li><a href="index.html">XMDS2 2.2.2 documentation</a> »</li>
</ul>
</div>
<div class="footer">
© Copyright 2008-2014, Graham Dennis, Joe Hope and Mattias Johnsson. Licensed under the GNU FDL.
Last updated on Oct 14, 2014.
Created using <a href="http://sphinx-doc.org/">Sphinx</a> 1.2b1.
</div>
</body>
</html>
|
